library(tidyverse)
library(plotrix)    #std.error
library(hypr)       #for generating contrast coding
library(ggthemes)
library(ggh4x)
library(dagitty)    #dagitty for causal diagrams
library(DiagrammeR) #grViz for causal diagrams
library(DiagrammeRsvg)
library(rsvg)
library(plotly)     #ggplotly for interactive plots
library(brms)       #bayesian regression
library(bayestestR) #p_direction, hdi
library(tidybayes)  #add_epred_draws
library(emmeans)
library(doParallel) #set the number of cores
library(ppcor)      #partial correlation

# # for power analysis
# library(lme4)
# library(mixedpower)


### Set global options
options(digits = 3) # set the default number of digits to 3
cores = as.integer(detectCores(logical = FALSE) * 0.7) # set the number of cores to use
registerDoParallel(cores=cores) # register the number of cores to use for parallel processing
options(mc.cores = cores)
options(brms.backend = "cmdstanr")  #this will speed up the model fitting

### MCMC options
niter = 20000  #number of iterations
nwu = 2000 #number of warmups

### Rmd settings
knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE)
### Custom functions

########### pp_check_each_round ############
pp_check_each_round <- function(m, data, round_info) {
  df_sub = filter(data, round == round_info)
  p = pp_check(m,
               type = "bars",
               ndraws = 100,
               newdata = df_sub) +
    coord_cartesian(xlim = c(0, 10)) +
    ggtitle(round_info)
  return(p)
}


########### pp_check_each_condition ############
pp_check_each_condition <- function(m, data, condition_info) {
  df_sub = filter(data, condition == condition_info)
  p = pp_check(m,
               type = "bars",
               ndraws = 100,
               newdata = df_sub) +
    coord_cartesian(xlim = c(0, 10)) +
    ggtitle(condition_info)
  return(p)
}


########### plot_posterior ############
plot_posterior <- function(model, model2=NA, model3=NA, 
                           interaction=FALSE, include_intercept=FALSE, 
                           xlim_cond=1.5, xlim_round=2, xlim_lex=0.15){
  ### extract the posterior draws
  posterior_beta1 <- model %>% 
    gather_draws(`b_.*`, regex = TRUE) %>% 
    mutate(intercept = str_detect(.variable, "Intercept"),
           component = ifelse(str_detect(.variable, ":"), "Interaction", 
                              ifelse(str_detect(.variable, "round"), "Round", 
                                     ifelse(str_detect(.variable, "Intercept"), "Intercept",
                                            ifelse(str_detect(.variable, "lex_align"), "N. lex align",
                                                   "Visibility")))))
  
  if (length(model2) == 1){ #if model2 is NA 
    posterior_beta = posterior_beta1
  } else {
    posterior_beta1 <- posterior_beta1 %>% 
      filter(.variable != "b_num_lex_align")
    posterior_beta2 <- model2 %>% 
      gather_draws(`b_.*`, regex = TRUE) %>% 
      filter(.variable == "b_num_lex_align") %>% 
      mutate(component = "N. lex align")
    posterior_beta3 <- model3 %>% 
      gather_draws(`b_.*`, regex = TRUE) %>% 
      filter(.variable == "b_condition_sumAO_Sym") %>% 
      mutate(component = "Visibility")
    posterior_beta = rbind(posterior_beta1, posterior_beta2, posterior_beta3)
  }
  
  posterior_beta = posterior_beta %>% 
    mutate(.variable = recode(.variable, 
                              "b_Intercept" = "Intercept",
                              "b_conditionAsymAV" = "SymAV--AsymAV",
                              "b_conditionAO" = "SymAV--AO",
                              "b_conditionAsym_Sym" = "AsymAV--SymAV",
                              "b_conditionAO_Asym" = "AO--AsymAV",
                              "b_condition_sumAO_Sym" = "AO--SymAV",
                              "b_num_lex_align" = "N. lex align",
                              "b_round_c" = "Centered round",
                              "b_log_round_c" = "Centered log(round)",
                              "b_roundR12" = "R1--R2",
                              "b_roundR23" = "R2--R3",
                              "b_roundR34" = "R3--R4",
                              "b_roundR45" = "R4--R5",
                              "b_roundR56" = "R5--R6",
                              "b_conditionAsymAV:round1" = "SymAV--AsymAV: R1--R2",
                              "b_conditionAsymAV:round2" = "SymAV--AsymAV: R2--R3",
                              "b_conditionAsymAV:round3" = "SymAV--AsymAV: R3--R4",
                              "b_conditionAsymAV:round4" = "SymAV--AsymAV: R4--R5",
                              "b_conditionAsymAV:round5" = "SymAV--AsymAV: R5--R6",
                              "b_conditionAO:round1" = "SymAV--AO: R1--R2",
                              "b_conditionAO:round2" = "SymAV--AO: R2--R3",
                              "b_conditionAO:round3" = "SymAV--AO: R3--R4",
                              "b_conditionAO:round4" = "SymAV--AO: R4--R5",
                              "b_conditionAO:round5" = "SymAV--AO: R5--R6",
                              "b_conditionAsym_Sym:round_c" = "Centered round: Asym--Sym",
                              "b_conditionAO_Sym:round_c" = "Centered round: AO--Sym",
                              "b_conditionAO_Asym:round_c" = "Centered round: AO--Asym",
                              "b_conditionAsym_Sym:log_round_c" = "Centered log(round): Asym--Sym",
                              "b_conditionAO_Sym:log_round_c" = "Centered log(round): AO--Sym",
                              "b_conditionAO_Asym:log_round_c" = "Centered log(round): AO--Asym"),
           .variable = factor(.variable,
                              levels = c("AO--SymAV", "AO--AsymAV", "AsymAV--SymAV", 
                                         "R1--R2", "R2--R3", "R3--R4", "R4--R5", "R5--R6",
                                         "N. lex align")),
           component = factor(component, 
                              levels = c("Intercept", "Visibility", "Round", 
                                         "Interaction", "N. lex align")))
  
  if (include_intercept == F){
    posterior_beta = posterior_beta %>% filter(component != "Intercept")
  }
  ### change variables if only main effects are plotted
  if (interaction == F) {
    posterior_beta = filter(posterior_beta, !str_detect(.variable, ":"))
    fill_manual_values = c("steelblue", "steelblue", "steelblue")
  } else{
    fill_manual_values = c("steelblue", "steelblue", "steelblue", "steelblue")
  }
  
  
  ### plot the posterior distributions
  p_posterior = ggplot(posterior_beta, 
                       aes(x = .value, y = fct_rev(.variable),
                           fill = component)) +
    geom_vline(xintercept = 0, size = 1) +
    stat_halfeye(aes(slab_alpha = intercept),
                 normalize = "panels", 
                 slab_alpha = .5,
                 .width = c(0.89, 0.95), 
                 point_interval = "median_hdi") +
    scale_fill_manual(values = fill_manual_values) +
    scale_slab_alpha_discrete(range = c(0.8, 0.4)) +
    guides(fill = "none", slab_alpha = "none") +
    labs(x = "Coefficient", y = "Effect") +
    theme_clean(base_size = 15) +
    theme(axis.text.x = element_text(colour = "black", size = 14),
          axis.text.y = element_text(colour = "black", size = 14),
          axis.title = element_text(size = 15, face = 'bold'),
          axis.title.x = element_text(vjust = -2),
          axis.title.y = element_text(vjust = 2),
          legend.position = "none",
          strip.text = element_text(size = 15, face = 'bold'),
          strip.background = element_blank(),
          panel.grid.major.x = element_line(color = "grey90", 
                                            linetype = "solid",
                                            size = 0.5),
          panel.grid.major.y = element_line(color = "grey90", 
                                            linetype = "solid",
                                            size = 0.5),
          plot.background = element_blank(),
          plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
    facet_wrap(vars(component), ncol = 3, scales = "free") + 
    facetted_pos_scales(
      x = list(
        scale_x_continuous(limits = c(-xlim_cond, xlim_cond)),
        scale_x_continuous(limits = c(-xlim_round, xlim_round)),
        scale_x_continuous(limits = c(-xlim_lex, xlim_lex),
                           breaks = c(-0.1, 0, 0.1))))
  
  return(p_posterior)
}


########### pp_update_plot ############
pp_update_plot <- function(post_sample, model_type="nb", interaction=TRUE){
  sum = ifelse("b_condition_sumAO_Sym" %in% colnames(post_sample), T, F)
  round_bd = ifelse("b_roundR12" %in% colnames(post_sample), T, F)
  
  intercept = ggplot(post_sample) +
    geom_density(aes(prior_Intercept), fill="steelblue", color="black",alpha=0.6) +
    geom_density(aes(b_Intercept), fill="#FC4E07", color="black",alpha=0.6) + 
    xlab('Intercept') +
    theme_classic()
  
  ### Visibility condition
  if (sum == F){
    cond1 = ggplot(post_sample) +
      geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
      geom_density(aes(b_conditionAsym_Sym), fill="#FC4E07", color="black",alpha=0.6) + 
      xlab('Asym--Sym') +
      theme_classic()
    cond2 = ggplot(post_sample) +
      geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
      geom_density(aes(b_conditionAO_Asym), fill="#FC4E07", color="black",alpha=0.6) + 
      xlab('AO--Asym') +
      theme_classic()
  } else {
    cond1 = ggplot(post_sample) +
      geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
      geom_density(aes(b_condition_sumAO_Sym), fill="#FC4E07", color="black",alpha=0.6) + 
      xlab('AO--Sym') +
      theme_classic()
    cond2 = ggplot(post_sample) +
      geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
      geom_density(aes(b_condition_sumAsym_Sym), fill="#FC4E07", color="black",alpha=0.6) + 
      xlab('Asym--Sym') +
      theme_classic()
  }
  
  ### Round
  if (interaction) {
    if (round_bd){
      r1 = ggplot(post_sample) +
        geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
        geom_density(aes(b_roundR12), fill="#FC4E07", color="black",alpha=0.6) + 
        xlab('R1--R2') +
        theme_classic()
      r2 = ggplot(post_sample) +
        geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
        geom_density(aes(b_roundR23), fill="#FC4E07", color="black",alpha=0.6) + 
        xlab('R2--R3') +
        theme_classic()
      r3 = ggplot(post_sample) +
        geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
        geom_density(aes(b_roundR34), fill="#FC4E07", color="black",alpha=0.6) + 
        xlab('R3--R4') +
        theme_classic()
      r4 = ggplot(post_sample) +
        geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
        geom_density(aes(b_roundR45), fill="#FC4E07", color="black",alpha=0.6) + 
        xlab('R4--R5') +
        theme_classic()
      r5 = ggplot(post_sample) +
        geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
        geom_density(aes(b_roundR56), fill="#FC4E07", color="black",alpha=0.6) + 
        xlab('R5--R6') +
        theme_classic()
    } else {
      if ("b_round_c" %in% colnames(post_sample)) {
        round = ggplot(post_sample) +
          geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
          geom_density(aes(b_round_c), fill="#FC4E07", color="black",alpha=0.6) + 
          xlab('Centered round') +
          theme_classic()
        if (sum == F){
          cond1_round = ggplot(post_sample) +
            geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
            geom_density(aes(`b_conditionAsym_Sym:round_c`), fill="#FC4E07", color="black",alpha=0.6) + 
            xlab('Centered Round: Asym--Sym') +
            theme_classic()
          cond2_round = ggplot(post_sample) +
            geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
            geom_density(aes(`b_conditionAO_Asym:round_c`), fill="#FC4E07", color="black",alpha=0.6) + 
            xlab('Centered Round: AO--Asym') +
            theme_classic()
        } else {
          cond1_round = ggplot(post_sample) +
            geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
            geom_density(aes(`b_condition_sumAO_Sym:round_c`), fill="#FC4E07", color="black",alpha=0.6) + 
            xlab('Centered Round: AO--Sym') +
            theme_classic()
          cond2_round = ggplot(post_sample) +
            geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
            geom_density(aes(`b_condition_sumAsym_Sym:round_c`), fill="#FC4E07", color="black",alpha=0.6) + 
            xlab('Centered Round: Asym--Sym') +
            theme_classic()
        }
      } else {
        round = ggplot(post_sample) +
          geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
          geom_density(aes(b_log_round_c), fill="#FC4E07", color="black",alpha=0.6) + 
          xlab('Centered log(round)') +
          theme_classic()
        if (sum == F){
          cond1_round = ggplot(post_sample) +
            geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
            geom_density(aes(`b_conditionAsym_Sym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) + 
            xlab('Centered log(round): Asym--Sym') +
            theme_classic()
          cond2_round = ggplot(post_sample) +
            geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
            geom_density(aes(`b_conditionAO_Asym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) + 
            xlab('Centered log(round): AO--Asym') +
            theme_classic()
        } else {
          cond1_round = ggplot(post_sample) +
            geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
            geom_density(aes(`b_condition_sumAO_Sym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) + 
            xlab('Centered log(round): AO--Sym') +
            theme_classic()
          cond2_round = ggplot(post_sample) +
            geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
            geom_density(aes(`b_condition_sumAsym_Sym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) + 
            xlab('Centered log(round): Asym--Sym') +
            theme_classic()
        }
      }
    }
  }
  
  if (model_type == "nb"){
    shape = ggplot(post_sample) +
      geom_density(aes(prior_shape), fill="steelblue", color="black",alpha=0.6) +
      geom_density(aes(shape), fill="#FC4E07", color="black",alpha=0.6) + 
      xlab('Shape') +
      scale_x_continuous(limits = c(0, 10)) +
      theme_classic()} 
  else if (model_type == "zinb") {
    shape = ggplot(post_sample) +
      geom_density(aes(prior_shape), fill="steelblue", color="black",alpha=0.6) +
      geom_density(aes(shape), fill="#FC4E07", color="black",alpha=0.6) + 
      xlab('Shape') +
      scale_x_continuous(limits = c(0, 10)) +
      theme_classic()
    zi = ggplot(post_sample) +
      geom_density(aes(prior_zi), fill="steelblue", color="black",alpha=0.6) +
      geom_density(aes(zi), fill="#FC4E07", color="black",alpha=0.6) + 
      xlab('Zero-inflation') +
      theme_classic()} 
  else if (model_type == "zibt"){
    phi = ggplot(post_sample) +
      geom_density(aes(prior_phi), fill="steelblue", color="black",alpha=0.6) +
      geom_density(aes(phi), fill="#FC4E07", color="black",alpha=0.6) + 
      xlab('Precision') +
      theme_classic()
    zoi = ggplot(post_sample) +
      geom_density(aes(prior_zoi), fill="steelblue", color="black",alpha=0.6) +
      geom_density(aes(zoi), fill="#FC4E07", color="black",alpha=0.6) + 
      xlab('Zero-inflation') +
      theme_classic()
    coi = ggplot(post_sample) +
      geom_density(aes(prior_coi), fill="steelblue", color="black",alpha=0.6) +
      geom_density(aes(coi), fill="#FC4E07", color="black",alpha=0.6) + 
      xlab('One-inflation') +
      theme_classic()
  }
  
  ### display the plots
  if (interaction==F){
    if (model_type=="nb"){
      gridExtra::grid.arrange(intercept, cond1, cond2, shape, ncol=2)} 
    else if (model_type=="zinb"){
      gridExtra::grid.arrange(intercept, cond1, cond2, shape, zi, ncol=3)} 
    else if (model_type=="zibt"){
      gridExtra::grid.arrange(intercept, cond1, cond2, phi, zoi, coi, ncol=3)
    }
  } else {
    if (round_bd == T){
      if (model_type=="nb"){
        gridExtra::grid.arrange(intercept, cond1, cond2, shape, ncol=3)} 
      else if (model_type=="zinb"){
        gridExtra::grid.arrange(intercept, cond1, cond2, shape, zi, ncol=3)} 
      else if (model_type=="zibt"){
        gridExtra::grid.arrange(intercept, cond1, cond2, phi, zoi, coi, ncol=3)
      }
    } else {
      if (model_type=="nb"){
        gridExtra::grid.arrange(intercept, cond1, cond2, round,
                                cond1_round, cond2_round, shape, ncol=3)} 
      else if (model_type=="zinb"){
        gridExtra::grid.arrange(intercept, cond1, cond2, round,
                                cond1_round, cond2_round, shape, zi, ncol=3)} 
      else if (model_type=="zibt"){
        gridExtra::grid.arrange(intercept, cond1, cond2, round,
                                cond1_round, cond2_round,  phi, zoi, coi, ncol=3)
      }
    }
  }
}

1 =====Data preparation=====

1.1 Load data

### trial_info.csv
df_trial_info = read.csv("data/trial_info.csv", stringsAsFactors = TRUE) %>% 
  filter(num_words != 0) %>%  # remove trials that were accidentally skipped) %>%
  mutate(pair = factor(pair),
         round_c = as.integer(round) - mean(as.integer(round)),
         round_n = factor(round),
         round = factor(paste0("R", round)),
         log_round = log(as.integer(round)),
         log_round_c = log_round - mean(log_round),
         condition = factor(condition,
                            levels = c("Sym", "Asym", "AO"),
                            labels = c("SymAV", "AsymAV", "AO")),
         condition_sum = condition,
         duration_s = duration_ms / 1000,
         n_words_100 = num_words / 100,
         n_iconic_per_100words = num_iconic_gestures / n_words_100,
         n_iconic_100 = num_iconic_gestures / 100) %>%
  dplyr::select(pair, condition, round, round_n, 
                round_c, log_round, log_round_c, 
                duration_ms, duration_s, duration_m,
                director:n_iconic_100) %>% 
  rename(trial_duration_s = duration_s, trial_duration_m = duration_m)

head(df_trial_info)
### condition info
df_condition_info = read.csv("data/condition_info.csv", stringsAsFactors = TRUE) %>% 
  mutate(pair = factor(pair),
         condition = factor(condition,
                            levels = c("Sym", "Asym", "AO"),
                            labels = c("SymAV", "AsymAV", "AO")))


2 Demographics

2.1 Participant

### participants
df_participant = read_tsv("data/participant_info.txt") %>% 
  mutate(pair = factor(as.integer(pair)),
         condition = factor(condition))

df_participant_used = df_participant %>% 
  filter(used == "Y") %>% 
  mutate(condition = factor(condition, levels = c("Sym", "Asym", "AO")))
df_participant_used
mean_age = mean(df_participant_used$age)
sd_age = sd(df_participant_used$age)
range_age = paste(range(df_participant_used$age)[1], range(df_participant_used$age)[2], sep = " - ")
n_participants = nrow(df_participant_used)
n_female = nrow(filter(df_participant_used, gender == "F"))
n_male = nrow(filter(df_participant_used, gender == "M"))

data.frame(mean_age, sd_age, range_age, n_participants, n_female, n_male)
### dyads
# make a dataframe with dyad information (for each pair, gender & age of speaker 1 & 2)
df_dyad = df_participant_used %>% 
  dplyr::select(-notes) %>%
  pivot_wider(names_from = speaker, 
              values_from = c("gender", "age")) %>% 
  mutate(mixed_gender = ifelse(gender_A == gender_B, "N", "Y"),
         mixed_gender = ifelse(is.na(mixed_gender), "Y", mixed_gender),
         gender_dyad = paste0(gender_A, "_", gender_B))

df_dyad %>% 
  group_by(condition, mixed_gender) %>% 
  summarize(n = n())
df_dyad %>% 
  group_by(condition, mixed_gender, gender_dyad) %>% 
  summarize(n = n())

2.2 Trial information & iconic gestures

summarize_demographics <- function(df) {
  df %>%  
    summarize(total_s = sum(trial_duration_s),
              total_m = total_s / 60,
              ### trial duration ###
              mean_trial_dur_s = mean(trial_duration_s),
              mean_trial_dur_m = mean(trial_duration_m),
              sd_trial_dur_m = sd(trial_duration_m),
              lci_trial_dur_m = mean_trial_dur_m - 1.96 * sd_trial_dur_m / sqrt(n()),
              uci_trial_dur_m = mean_trial_dur_m + 1.96 * sd_trial_dur_m / sqrt(n()),
              ### words ###
              # number of words
              num_words_total = sum(num_words),
              mean_num_words = mean(num_words),
              num_words_100 = mean(num_words / 100),
              sd_num_words = sd(num_words),
              lci_num_words = mean_num_words - 1.96 * sd_num_words / sqrt(n()),
              uci_num_words = mean_num_words + 1.96 * sd_num_words / sqrt(n()),
              num_words_per_min = num_words_total / total_m,
              # number of content words
              num_content_total = sum(num_content_words),
              mean_num_content = mean(num_content_words),
              sd_num_content = sd(num_content_words),
              lci_num_content = mean_num_content - 1.96 * sd_num_content / sqrt(n()),
              uci_num_content = mean_num_content + 1.96 * sd_num_content / sqrt(n()),
              num_content_per_min = num_content_total / total_m,
              ### lexical alignment ###
              # raw frequency
              num_lex_align_total = sum(num_lex_align, na.rm = T),
              mean_num_lex_align = mean(num_lex_align, na.rm = T),
              sd_num_lex_align = sd(num_lex_align, na.rm = T),
              lci_num_lex_align = mean_num_lex_align - 1.96 * sd_num_lex_align / sqrt(n()),
              uci_num_lex_align = mean_num_lex_align + 1.96 * sd_num_lex_align / sqrt(n()),
              num_lex_align_per_min = num_lex_align_total / total_m,
              # rate per 100 words
              mean_lex_align_per_100words = mean(num_lex_align / num_words_100, na.rm=T),
              sd_lex_align_per_100words = sd(num_lex_align / num_words_100, na.rm=T),
              se_lex_align_per_100words = std.error(num_lex_align / num_words_100, na.rm=T),
              lci_lex_align_per_100words = mean_lex_align_per_100words - 1.96 * se_lex_align_per_100words,
              uci_lex_align_per_100words = mean_lex_align_per_100words + 1.96 * se_lex_align_per_100words,
              ### iconic gestures ###
              # raw frequency
              num_iconic_total = sum(num_iconic_gestures, na.rm = T),
              mean_num_iconic = mean(num_iconic_gestures, na.rm = T),
              sd_num_iconic = sd(num_iconic_gestures, na.rm = T),
              lci_num_iconic = mean_num_iconic - 1.96 * sd_num_iconic / sqrt(n()),
              uci_num_iconic = mean_num_iconic + 1.96 * sd_num_iconic / sqrt(n()),
              num_iconic_per_min = num_iconic_total / total_m,
              # rate per 100 words
              mean_iconic_per_100words = mean(num_iconic_gestures / num_words_100, na.rm=T),
              sd_iconic_per_100words = sd(num_iconic_gestures / num_words_100, na.rm=T),
              se_iconic_per_100words = std.error(num_iconic_gestures / num_words_100, na.rm=T),
              lci_iconic_per_100words = mean_iconic_per_100words - 1.96 * se_iconic_per_100words,
              uci_iconic_per_100words = mean_iconic_per_100words + 1.96 * se_iconic_per_100words,
              ### gestural alignment ###
              # raw frequency
              num_gest_align_total = sum(num_gestural_align, na.rm = T),
              mean_num_gest_align = mean(num_gestural_align, na.rm = T),
              sd_num_gest_align = sd(num_gestural_align, na.rm = T),
              lci_num_gest_align = mean_num_gest_align - 1.96 * sd_num_gest_align / sqrt(n()),
              uci_num_gest_align = mean_num_gest_align + 1.96 * sd_num_gest_align / sqrt(n()),
              num_gest_align_per_min = num_gest_align_total / total_m,
              # rate per 100 words
              mean_gest_align_per_100words = mean(num_gestural_align / num_words_100, na.rm=T),
              sd_gest_align_per_100words = sd(num_gestural_align / num_words_100, na.rm=T),
              se_gest_align_per_100words = std.error(num_gestural_align / num_words_100, na.rm=T),
              lci_gest_align_per_100words = mean_gest_align_per_100words - 1.96 * se_gest_align_per_100words,
              uci_gest_align_per_100words = mean_gest_align_per_100words + 1.96 * se_gest_align_per_100words,
              # rate per iconic gestures
              mean_gest_align_prop = mean(num_gestural_align / num_iconic_gestures, na.rm=T),
              sd_gest_align_per_prop = sd(num_gestural_align / num_iconic_gestures, na.rm=T),
              se_gest_align_per_prop = std.error(num_gestural_align / num_iconic_gestures, na.rm=T),
              lci_gest_align_per_prop = mean_gest_align_prop - 1.96 * se_gest_align_per_prop,
              uci_gest_align_per_prop = mean_gest_align_prop + 1.96 * se_gest_align_per_prop,
              ### number of trials ###
              trial_n = n()) %>% 
    ungroup()
}

summarize_dur <- function(df){
  df %>% 
    summarize(mean_dur_m = mean(total_m),
              sd_dur_m = sd(total_m),
              se_dur_m = std.error(total_m),
              lci_dur_m = mean_dur_m - 1.96 * se_dur_m,
              uci_dur_m = mean_dur_m + 1.96 * se_dur_m) %>% 
    ungroup()
}


2.2.1 mean by condition

### demographics by pair
trial_info_pair = df_trial_info %>% 
  group_by(pair) %>% 
  summarize_demographics() %>% 
  left_join(., df_condition_info) %>% 
  dplyr::select(pair, condition, total_s:trial_n)

### calculate mean experiment duration
mean_dur_cond = trial_info_pair %>% 
  group_by(condition) %>% 
  summarize_dur()

### summary statistics
trial_info_cond  = df_trial_info %>% 
  group_by(condition) %>% 
  summarize_demographics() %>% 
  left_join(., mean_dur_cond) %>% 
  dplyr::select(condition, total_m, mean_dur_m:uci_dur_m, 
                everything(), 
                -ends_with("_s"))

trial_info_cond


2.2.2 mean by round

trial_info_pair_round = df_trial_info %>% 
  group_by(pair, round, round_n) %>%
  summarize_demographics() %>% 
  left_join(., df_condition_info)

### calculate mean round duration
mean_dur_round = trial_info_pair_round %>% 
  group_by(round, round_n) %>% 
  summarize_dur()

trial_info_round = df_trial_info %>% 
  group_by(round, round_n) %>% 
  summarize_demographics() %>% 
  left_join(., mean_dur_round) %>%
  dplyr::select(round, round_n, total_m, mean_dur_m:uci_dur_m, 
                everything(), 
                -ends_with("_s"))

trial_info_round


2.2.3 mean by condition x round

### calculate mean duration by condition x round
mean_dur_cond_round = trial_info_pair_round %>% 
  group_by(condition, round, round_n) %>% 
  summarize_dur()

trial_info_cond_round = df_trial_info %>% 
  group_by(condition, round, round_n) %>% 
  summarize_demographics() %>% 
  left_join(., mean_dur_cond_round) %>%
  dplyr::select(condition, round, round_n,
                total_m, mean_dur_m:uci_dur_m, 
                everything(), 
                -ends_with("_s"))

trial_info_cond_round



3 =====Contrast coding=====

### visibility condition: difference coding
h_cond = hypr(AO_Asym = AsymAV ~ AO,
              Asym_Sym = SymAV ~ AsymAV,
              levels = levels(df_trial_info$condition))
h_cond
## hypr object containing 2 null hypotheses:
##  H0.AO_Asym: 0 = AsymAV - AO
## H0.Asym_Sym: 0 = SymAV - AsymAV
## 
## Call:
## hypr(AO_Asym = ~AsymAV - AO, Asym_Sym = ~SymAV - AsymAV, levels = c("SymAV", 
## "AsymAV", "AO"))
## 
## Hypothesis matrix (transposed):
##        AO_Asym Asym_Sym
## SymAV   0       1      
## AsymAV  1      -1      
## AO     -1       0      
## 
## Contrast matrix:
##        AO_Asym Asym_Sym
## SymAV   1/3     2/3    
## AsymAV  1/3    -1/3    
## AO     -2/3    -1/3
contrasts(df_trial_info$condition) = contr.hypothesis(h_cond)


### visibility condition: sum coding
h_cond = hypr(AO_Sym = SymAV ~ AO,
              Asym_Sym = SymAV ~ AsymAV,
              levels = levels(df_trial_info$condition))
h_cond
## hypr object containing 2 null hypotheses:
##   H0.AO_Sym: 0 = SymAV - AO
## H0.Asym_Sym: 0 = SymAV - AsymAV
## 
## Call:
## hypr(AO_Sym = ~SymAV - AO, Asym_Sym = ~SymAV - AsymAV, levels = c("SymAV", 
## "AsymAV", "AO"))
## 
## Hypothesis matrix (transposed):
##        AO_Sym Asym_Sym
## SymAV   1      1      
## AsymAV  0     -1      
## AO     -1      0      
## 
## Contrast matrix:
##        AO_Sym Asym_Sym
## SymAV   1/3    1/3    
## AsymAV  1/3   -2/3    
## AO     -2/3    1/3
contrasts(df_trial_info$condition_sum) = contr.hypothesis(h_cond)


### round
bacward_diff = hypr(R12 = R2 ~ R1,
                    R23 = R3 ~ R2,
                    R34 = R4 ~ R3,
                    R45 = R5 ~ R4,
                    R56 = R6 ~ R5,
                    levels = levels(df_trial_info$round))
bacward_diff
## hypr object containing 5 null hypotheses:
## H0.R12: 0 = R2 - R1
## H0.R23: 0 = R3 - R2
## H0.R34: 0 = R4 - R3
## H0.R45: 0 = R5 - R4
## H0.R56: 0 = R6 - R5
## 
## Call:
## hypr(R12 = ~R2 - R1, R23 = ~R3 - R2, R34 = ~R4 - R3, R45 = ~R5 - 
##     R4, R56 = ~R6 - R5, levels = c("R1", "R2", "R3", "R4", "R5", 
## "R6"))
## 
## Hypothesis matrix (transposed):
##    R12 R23 R34 R45 R56
## R1 -1   0   0   0   0 
## R2  1  -1   0   0   0 
## R3  0   1  -1   0   0 
## R4  0   0   1  -1   0 
## R5  0   0   0   1  -1 
## R6  0   0   0   0   1 
## 
## Contrast matrix:
##    R12  R23  R34  R45  R56 
## R1 -5/6 -2/3 -1/2 -1/3 -1/6
## R2  1/6 -2/3 -1/2 -1/3 -1/6
## R3  1/6  1/3 -1/2 -1/3 -1/6
## R4  1/6  1/3  1/2 -1/3 -1/6
## R5  1/6  1/3  1/2  2/3 -1/6
## R6  1/6  1/3  1/2  2/3  5/6
contrasts(df_trial_info$round) = contr.hypothesis(bacward_diff)

4 =====Number of iconic gestures=====

4.1 —DataViz: frequency—

4.1.1 bp: total by condition

bp_iconic_by_cond = ggplot(data=trial_info_pair, 
                           aes(x=condition, y=num_iconic_total, fill=condition)) +
  geom_jitter(aes(color = pair), 
              size = 1, alpha = 1, 
              width = 0.07, height = 0) +
  geom_boxplot(width = .2,
               outlier.shape = NA, alpha = 0.7) +
  scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
  labs(x = "Visibility",
       y = "Total N of iconic gestures per pair") +
  theme_clean(base_size = 15) +
  theme(axis.text.x = element_text(colour = "black", size = 14),
        axis.text.y = element_text(colour = "black", size = 14),
        axis.title = element_text(size = 15, face = 'bold'),
        axis.title.x = element_text(vjust = -2),
        axis.title.y = element_text(vjust = 2),
        legend.position = "none",
        strip.text = element_text(size = 15, face = 'bold'),
        plot.background = element_blank(),
        plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))

ggplotly(bp_iconic_by_cond)


4.1.2 bp: mean by condition

bp_mean_iconic_by_cond = ggplot(data=trial_info_pair, 
                                aes(x=condition, y=mean_num_iconic, fill=condition)) +
  geom_jitter(aes(color = pair), 
              size = 1, alpha = 1, 
              width = 0.07, height = 0) +
  geom_boxplot(width = .2,
               outlier.shape = NA, alpha = 0.7) +
  geom_point(data = trial_info_cond, 
             aes(y = mean_num_iconic), 
             size = 2, shape = 4) +
  scale_y_continuous(limits = c(0, 4), breaks = seq(0, 4, 1)) +
  scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
  labs(x = "Visibility",
       y = "Mean N of iconic gestures per trial") +
  theme_clean(base_size = 15) +
  theme(axis.text.x = element_text(colour = "black", size = 14),
        axis.text.y = element_text(colour = "black", size = 14),
        axis.title = element_text(size = 15, face = 'bold'),
        axis.title.x = element_text(vjust = -2),
        axis.title.y = element_text(vjust = 2),
        legend.position = "none",
        strip.text = element_text(size = 15, face = 'bold'),
        plot.background = element_blank(),
        plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))

ggplotly(bp_mean_iconic_by_cond)


4.1.3 bp: mean by condition x round

bp_mean_iconic_by_round_cond = ggplot(data=df_trial_info, 
                                      aes(x=round, y=num_iconic_gestures, fill=condition)) +
  geom_boxplot(outlier.shape = NA,
               alpha = 0.7) +
  scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
  scale_y_continuous(limits = c(0, 15), breaks = seq(0, 15, 5)) +
  labs(x = "Round",
       y = "Mean N of iconic gestures per trial") +
  theme_clean(base_size = 15) +
  theme(axis.text.x = element_text(colour = "black", size = 14),
        axis.text.y = element_text(colour = "black", size = 14),
        axis.title = element_text(size = 15, face = 'bold'),
        axis.title.x = element_text(vjust = -2),
        axis.title.y = element_text(vjust = 2),
        legend.position = "top",
        strip.text = element_text(size = 15, face = 'bold'),
        plot.background = element_blank(),
        plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))

bp_mean_iconic_by_round_cond

The figure shows that the decline in iconic gestures across rounds does not follow a linear pattern. This suggests that log-transformed round may improve the model fit.


4.2 —Zero-inflated Poisson regression models—

4.2.1 Prior specification

We will use weakly informative priors for the regression each parameter. To specify priors for the intercept, which represents the expected number of iconic gestures in the SymAV condition, we will refer to the number of iconic gestures reported in Akamine et al. (2024). In the paper, the authors analyzed data from 19 dyads involving in the same task as the current study but in co-present interaction. The total number of iconic gestures reported was 4413, which was collected from 19 dyads, each performing 96 trials. Therefore, the expected number of iconic gestures per trial is \(4413 / (19 * 96) = 2.42\).

As Poisson regression models use a log-link, we need to specify priors on the log scale. Taking the log of 2.42, we get 0.88. In order to allow data to inform the posterior, we will use a normal prior with a mean of 0.88 and a relatively wide standard deviation of 3.

For the fixed effects, we will set unbiased priors with a mean of 0 and a standard deviation of 2.

For the standard deviation of the random effects, we set the prior to be normal with mean 0 and standard deviation 2. For the correlation between the random effects, we set the prior to be LKJ(2).

For the models including the round as fixed effects (whether backward-difference coded or centered + log-transformed), the intercept will represent the mean expected number of iconic gestures (ground mean) in the SymAV condition. As the meaning of the intercept doesn’t change when adding the round variable, we use the same prior for the intercept.

priors_rint = c(
  prior(normal(0.88, 2), class = Intercept),
  prior(normal(0, 2), class = b),
  prior(normal(0, 2), class = sd))

priors_rslope = c(
  prior(normal(0.88, 2), class = Intercept),
  prior(normal(0, 2), class = b),
  prior(normal(0, 2), class = sd),
  prior(lkj(2), class = cor))


4.2.2 Model 1: condition (pair)

First, we will model the number of iconic gestures per trial with fixed effects for condition. As we expect no gestures in many trials based on previous study (Akamine et al., 2024), we will use a zero-inflated Poisson regression model with a log link function.


4.2.2.1 Prior predictive check

mp_iconic_cond_prior = brm(num_iconic_gestures ~ 1 + condition + 
                             (1|pair) + (1|target),
                           family = poisson(),
                           prior = priors_rint,
                           sample_prior = "only",
                           data = df_trial_info,
                           file = "models/mp_iconic_cond_prior")

pp_check(mp_iconic_cond_prior, ndraws = 100, type = "bars") +
  coord_cartesian(xlim = c(0, 20),
                  ylim = c(0, 4000))

The prior predictive check shows that the model generates data that are similar to the observed data. The model generates a far wider range of data than the observed data, allowing the model to be informed by the data.

Next, we will fit the model to the observed data.


4.2.2.2 Fit the model

mp_iconic_cond = brm(num_iconic_gestures ~ 1 + condition + 
                       (1|pair) + (1|target),
                     family = zero_inflated_poisson(),
                     prior = priors_rint,
                     data = df_trial_info,
                     sample_prior = T,
                     warmup = nwu, iter = niter,
                     control = list(adapt_delta = 0.9, 
                                    max_treedepth = 15),
                     file = "models/mp_iconic_cond")

summary(mp_iconic_cond)
##  Family: zero_inflated_poisson 
##   Links: mu = log; zi = identity 
## Formula: num_iconic_gestures ~ 1 + condition + (1 | pair) + (1 | target) 
##    Data: df_trial_info (Number of observations: 4316) 
##   Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
##          total post-warmup draws = 72000
## 
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.55      0.07     0.44     0.69 1.00    10772    22264
## 
## ~target (Number of levels: 16) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.21      0.05     0.14     0.33 1.00    12904    27375
## 
## Regression Coefficients:
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             0.91      0.10     0.71     1.11 1.00     8817    16940
## conditionAO_Asym      0.29      0.21    -0.12     0.70 1.00     8349    15439
## conditionAsym_Sym    -0.07      0.20    -0.48     0.33 1.00     8526    15481
## 
## Further Distributional Parameters:
##    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## zi     0.41      0.01     0.39     0.43 1.00    56110    51525
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayestestR::hdi(mp_iconic_cond)
bayestestR::hdi(mp_iconic_cond, ci = 0.89)


4.2.2.3 Posterior predictive check

pp_check_sym = pp_check_each_condition(mp_iconic_cond, df_trial_info, "SymAV")
pp_check_asym = pp_check_each_condition(mp_iconic_cond, df_trial_info, "AsymAV")
pp_check_ao = pp_check_each_condition(mp_iconic_cond, df_trial_info, "AO")

gridExtra::grid.arrange(pp_check_sym, pp_check_asym, pp_check_ao, ncol = 2)

The posterior predictive check shows that the model generates data that are similar to the observed data. This suggests that the model is a good fit to the data.


4.2.2.4 Model Effect

Here, we will compare which regression model is suitable for the data, poisson or negative binomial. We will follow the procedure introduced in Winter & Burkner (2021).

First thing first, we will fit the negative binomial regression model to the data. Then we compare the predictive power of the models using loo_compare() function.

mnb_iconic_cond = brm(num_iconic_gestures ~ 1 + condition + 
                        (1|pair) + (1|target),
                      family = negbinomial(),
                      prior = priors_rint,
                      data = df_trial_info,
                      sample_prior = T,
                      warmup = nwu, iter = niter,
                      control = list(adapt_delta = 0.9, 
                                     max_treedepth = 15),
                      file = "models/mnb_iconic_cond")

mnb_iconic_cond_zi = brm(num_iconic_gestures ~ 1 + condition + 
                           (1|pair) + (1|target),
                         family = zero_inflated_negbinomial(),
                         prior = priors_rint,
                         data = df_trial_info,
                         sample_prior = T,
                         warmup = nwu, iter = niter,
                         control = list(adapt_delta = 0.9, 
                                        max_treedepth = 15),
                         file = "models/mnb_iconic_cond_zi")

### loo compare
if (!file.exists("models/loo_pois_iconic.rds")){
  pois_loo = loo(mp_iconic_cond)
  saveRDS(pois_loo, file = "models/loo_pois_iconic.rds")
}

if (!file.exists("models/loo_nb_iconic.rds")){
  nb_loo = loo(mnb_iconic_cond)
  saveRDS(nb_loo, file = "models/loo_nb_iconic.rds")
}

if (!file.exists("models/loo_zinb_iconic.rds")){
  zinb_loo = loo(mnb_iconic_cond_zi)
  saveRDS(zinb_loo, file = "models/loo_zinb_iconic.rds")
}

pois_loo = readRDS("models/loo_pois_iconic.rds")
nb_loo = readRDS("models/loo_nb_iconic.rds")
zinb_loo = readRDS("models/loo_zinb_iconic.rds")

loo_compare(pois_loo, nb_loo, zinb_loo)
##                    elpd_diff se_diff
## mnb_iconic_cond       0.0       0.0 
## mnb_iconic_cond_zi   -0.9       1.5 
## mp_iconic_cond     -947.7     102.0

The loo_compare() function shows that the negative binomial regression has a larger predictive power (elpd_diff) and smaller standard error (se_diff) than the zero-inflated negative binomial regression and the Poisson regression. Therefore, we will use the negative binomial regression model for further analyses.



4.3 —Negative binomial regression models—

4.3.1 Model 2: condition x round (pair)

Next, we will run the negative binomial regression model to predict the number of iconic gestures per trial by condition and round. We will use the priors_round for the model.


4.3.1.1 Model effect (log round)

We saw earlier that the decline in # iconic gestures doesn’t follow a linear pattern: The decline is stepper in the initial rounds and then plateaus afterwards. To capture this trend better without making 5 levels for Round (which is the case for bd coding), we will check if the model performance is comparable for backward-difference coded round and centered log-transformed round.

mnb_iconic_cond_round = brm(num_iconic_gestures ~ 1 + condition * round + 
                              (1+round|pair) + (1|target),
                            family = negbinomial(),
                            prior = priors_rslope,
                            data = df_trial_info,
                            sample_prior = T,
                            warmup = nwu, iter = niter,
                            control = list(adapt_delta = 0.9, 
                                           max_treedepth = 15),
                            file = "models/mnb_iconic_cond_round")

mnb_iconic_cond_round_c = brm(num_iconic_gestures ~ 1 + condition * round_c + 
                                (1+round_c|pair) + (1|target),
                              family = negbinomial(),
                              prior = priors_rslope,
                              data = df_trial_info,
                              sample_prior = T,
                              warmup = nwu, iter = niter,
                              control = list(adapt_delta = 0.9, 
                                             max_treedepth = 15),
                              file = "models/mnb_iconic_cond_round_c")

mnb_iconic_cond_round_log = brm(num_iconic_gestures ~ 1 + condition * log_round_c + 
                                  (1+log_round_c|pair) + (1|target),
                                family = negbinomial(),
                                prior = priors_rslope,
                                data = df_trial_info,
                                sample_prior = T,
                                warmup = nwu, iter = niter,
                                control = list(adapt_delta = 0.9, 
                                               max_treedepth = 15),
                                file = "models/mnb_iconic_cond_round_log")


### loo compare
if (!file.exists("models/loo_nb_iconic_cond_round.rds")){
  nb_cond_round_loo = loo(mnb_iconic_cond_round)
  saveRDS(nb_cond_round_loo, file = "models/loo_nb_iconic_cond_round.rds")
}

if (!file.exists("models/loo_nb_iconic_cond_round_c.rds")){
  nb_cond_round_c_loo = loo(mnb_iconic_cond_round_c)
  saveRDS(nb_cond_round_c_loo, file = "models/loo_nb_iconic_cond_round_c.rds")
}

if (!file.exists("models/loo_nb_iconic_cond_round_log.rds")){
  nb_cond_round_log_loo = loo(mnb_iconic_cond_round_log)
  saveRDS(nb_cond_round_log_loo, file = "models/loo_nb_iconic_cond_round_log.rds")
}

nb_cond_round_loo = readRDS("models/loo_nb_iconic_cond_round.rds")
nb_cond_round_c_loo = readRDS("models/loo_nb_iconic_cond_round_c.rds")
nb_cond_round_log_loo = readRDS("models/loo_nb_iconic_cond_round_log.rds")

loo_compare(nb_cond_round_loo, nb_cond_round_c_loo, nb_cond_round_log_loo)
##                           elpd_diff se_diff
## mnb_iconic_cond_round_log   0.0       0.0  
## mnb_iconic_cond_round     -13.9       5.3  
## mnb_iconic_cond_round_c   -38.2      12.8

The leave-one-out (LOO) Effect shows that centered log-transformed round provides a larger predictive power (elpd_diff) and smaller standard error (se_diff) compared to the backward-difference coded round or centered round. This is in line with our observation that the decrease in the number of gestures does not follow a linear pattern, which explains why the centered round is disfavored here. Therefore, we will use the centered log round.


4.3.1.2 Prior predictive check

mnb_iconic_cond_round_log = brm(num_iconic_gestures ~ 1 + condition * log_round_c + 
                                  (1+log_round_c|pair) + (1|target),
                                family = negbinomial(),
                                prior = priors_rslope,
                                sample_prior = "only",
                                data = df_trial_info,
                                file = "models/mnb_iconic_cond_round_log_prior")

pp_check(mnb_iconic_cond_round_log, ndraws = 100, type = "bars") +
  coord_cartesian(xlim = c(0, 20),
                  ylim = c(0, 4000))

The prior predictive check shows that the model generates data that are similar to the observed data.


4.3.1.3 Fit the model

mnb_iconic_cond_round_log = brm(num_iconic_gestures ~ 1 + condition * log_round_c + 
                                  (1+log_round_c|pair) + (1|target),
                                family = negbinomial(),
                                prior = priors_round,
                                data = df_trial_info,
                                sample_prior = T,
                                warmup = nwu, iter = niter,
                                control = list(adapt_delta = 0.9, 
                                               max_treedepth = 15),
                                file = "models/mnb_iconic_cond_round_log")

summary(mnb_iconic_cond_round_log)
##  Family: negbinomial 
##   Links: mu = log; shape = identity 
## Formula: num_iconic_gestures ~ 1 + condition * log_round_c + (1 + log_round_c | pair) + (1 | target) 
##    Data: df_trial_info (Number of observations: 4316) 
##   Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
##          total post-warmup draws = 72000
## 
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.85      0.10     0.68     1.07 1.00    15511
## sd(log_round_c)                0.41      0.06     0.31     0.54 1.00    24199
## cor(Intercept,log_round_c)     0.64      0.11     0.38     0.82 1.00    30037
##                            Tail_ESS
## sd(Intercept)                 26870
## sd(log_round_c)               37632
## cor(Intercept,log_round_c)    44753
## 
## ~target (Number of levels: 16) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.30      0.07     0.20     0.45 1.00    18463    33962
## 
## Regression Coefficients:
##                               Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                        -0.05      0.15    -0.34     0.25 1.00
## conditionAO_Asym                  0.37      0.31    -0.24     0.98 1.00
## conditionAsym_Sym                -0.16      0.31    -0.77     0.45 1.00
## log_round_c                      -1.37      0.07    -1.51    -1.23 1.00
## conditionAO_Asym:log_round_c      0.00      0.17    -0.33     0.34 1.00
## conditionAsym_Sym:log_round_c    -0.10      0.17    -0.43     0.23 1.00
##                               Bulk_ESS Tail_ESS
## Intercept                        12579    22655
## conditionAO_Asym                 11906    21005
## conditionAsym_Sym                11866    21349
## log_round_c                      18656    32465
## conditionAO_Asym:log_round_c     19212    31521
## conditionAsym_Sym:log_round_c    19007    32931
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     2.06      0.13     1.83     2.32 1.00    81667    50845
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayestestR::hdi(mnb_iconic_cond_round_log)
# bayestestR::hdi(mnb_iconic_cond_round_log, ci = 0.89)

Based on the model estimates, the condition does not have a significant effect on the number of iconic gestures produced per trial. However, the round has a significant effect on the number of iconic gestures produced per trial: the number of iconic gestures significantly decreases as the round progresses.


4.3.1.4 Visualize the posterior distributions

# p = plot_posterior(mnb_iconic_cond_round_log, interaction=T, xlim=1)
# p


4.3.1.5 Probability of direction

p_direction(mnb_iconic_cond_round_log)

The probability of direction demonstrated significant effects for the round only.


4.3.1.6 Prior-posterior update plot

### condition
#Sample the parameters of interest:
posterior <- as_draws_df(mnb_iconic_cond_round_log)

#Plot the prior-posterior update plot for the estimated parameters:
pp_update_intercept = ggplot(posterior) +
  geom_density(aes(prior_Intercept), fill="steelblue", color="black",alpha=0.6) +
  geom_density(aes(b_Intercept), fill="#FC4E07", color="black",alpha=0.6) + 
  xlab('Intercept') +
  theme_classic()

pp_update_b1 = ggplot(posterior) +
  geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
  geom_density(aes(b_conditionAO_Asym), fill="#FC4E07", color="black",alpha=0.6) + 
  xlab('AO --Asym') +
  theme_classic()

pp_update_b2 = ggplot(posterior) +
  geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
  geom_density(aes(b_conditionAsym_Sym), fill="#FC4E07", color="black",alpha=0.6) + 
  xlab('Asym--Sym') +
  theme_classic()


### round
pp_update_round = ggplot(posterior) +
  geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
  geom_density(aes(b_log_round_c), fill="#FC4E07", color="black",alpha=0.6) + 
  xlab('Centered log(round)') +
  theme_classic()


gridExtra::grid.arrange(pp_update_intercept, pp_update_b1, pp_update_b2, pp_update_round)


4.3.1.7 Posterior predictive check

pp_check_sym = pp_check_each_condition(mnb_iconic_cond_round_log, df_trial_info, "SymAV")
pp_check_asym = pp_check_each_condition(mnb_iconic_cond_round_log, df_trial_info, "AsymAV")
pp_check_ao = pp_check_each_condition(mnb_iconic_cond_round_log, df_trial_info, "AO")

gridExtra::grid.arrange(pp_check_sym, pp_check_asym, pp_check_ao, ncol = 2)

pp_check_r1 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R1")
pp_check_r2 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R2")
pp_check_r3 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R3")
pp_check_r4 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R4")
pp_check_r5 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R5")
pp_check_r6 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R6")

gridExtra::grid.arrange(pp_check_r1, pp_check_r2, pp_check_r3, pp_check_r4, pp_check_r5, pp_check_r6, ncol = 3)


4.3.1.8 Pair-wise Effect

emmeans(mnb_iconic_cond_round_log, pairwise ~ condition)$contrasts
##  contrast       estimate lower.HPD upper.HPD
##  SymAV - AsymAV   -0.160    -0.769     0.454
##  SymAV - AO        0.214    -0.402     0.830
##  AsymAV - AO       0.375    -0.240     0.981
## 
## Point estimate displayed: median 
## Results are given on the log (not the response) scale. 
## HPD interval probability: 0.95
# emmeans(mnb_iconic_cond_round_log, pairwise ~ condition, level = 0.89)$contrasts


4.3.1.9 Visualize the model estimates

plot(conditional_effects(mnb_iconic_cond_round_log), ask = FALSE)



5 =====<Chosen> Rate of iconic gestures per 100 words=====

5.1 —DataViz: rate per 100 words—

5.1.1 bp: mean by condition

bp_mean_iconic_rate_by_cond = ggplot(data=trial_info_pair, 
                                     aes(x=condition, y=mean_iconic_per_100words, fill=condition)) +
  geom_jitter(aes(color = pair), 
              size = 1, alpha = 1, 
              width = 0.07, height = 0) +
  geom_boxplot(width = .2,
               outlier.shape = NA, alpha = 0.7) +
  geom_point(data = trial_info_cond, 
             aes(y = mean_iconic_per_100words), 
             shape = 21, size = 3, fill = "white") +
  scale_y_continuous(limits = c(0, 11), breaks = seq(0, 10, 2)) +
  scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
  labs(x = "Visibility",
       y = "Mean rate of iconic gestures") +
  theme_clean(base_size = 15) +
  theme(axis.text.x = element_text(colour = "black", size = 14),
        axis.text.y = element_text(colour = "black", size = 14),
        axis.title = element_text(size = 15, face = 'bold'),
        axis.title.x = element_text(vjust = -2),
        axis.title.y = element_text(vjust = 2),
        legend.position = "none",
        strip.text = element_text(size = 15, face = 'bold'),
        plot.background = element_blank(),
        plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))

ggplotly(bp_mean_iconic_rate_by_cond)


5.1.2 bp: mean by condition x round

pd = position_dodge(width = .75)

ggplot(data=df_trial_info, 
       aes(x=round, y=n_iconic_per_100words, fill=condition)) +
  geom_boxplot(outlier.shape = NA,
               alpha = 0.7) +
  geom_point(data = trial_info_cond_round, 
             aes(y = mean_iconic_per_100words, 
                 group = condition),
             position = pd,
             shape = 21, size = 2, fill = "white") +
  scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
  scale_y_continuous(limits = c(0, 22), breaks = seq(0, 20, 5)) +
  labs(x = "Round",
       y = "Mean rate of iconic gestures") +
  theme_clean(base_size = 15) +
  theme(axis.text.x = element_text(colour = "black", size = 14),
        axis.text.y = element_text(colour = "black", size = 14),
        axis.title = element_text(size = 15, face = 'bold'),
        axis.title.x = element_text(vjust = -2),
        axis.title.y = element_text(vjust = 2),
        legend.position = "top",
        strip.text = element_text(size = 15, face = 'bold'),
        plot.background = element_blank(),
        plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))



5.2 —<Chosen> ZI negative binomial regression models—

5.2.1 Prior specification

As the unit of the dependent variable is different from the previous model, we will set different priors for the rate of iconic gestures per 100 words. We will use weakly informative priors for the regression each parameter. To specify priors for the intercept, we will refer to the rate of iconic gestures per 100 words reported in Akamine et al. (2024). In the paper, the authors analyzed data from 19 dyads involving in the same task as the current study but in co-present interaction. The total number of iconic gestures reported was 4413, which was collected from 19 dyads, each performing 96 trials. The total number of words produced was 71695 (28152 content) words. Therefore, the expected rate of iconic gestures per 100 words is \(4413 / (71695 / 100) = 6.16\) (log(6.16) = 1.82).

For the fixed effects, we will set unbiased priors with a mean of 0 and a standard deviation of 0.5.

For the standard deviation of the random effects, we set the prior to be normal with mean 0 and standard deviation 0.5. For the correlation between the random effects, we set the prior to be LKJ(2).

For the models including the round as fixed effects (whether backward-difference coded or centered + log-transformed), the intercept will represent the mean expected number of iconic gestures (ground mean). As the meaning of the intercept doesn’t change when adding the round variable, we use the same prior for the intercept.

Note that we do not expect the rate of iconic gestures to change across rounds (i.e., we expect the number of words and gestures to decrease at an approximately same pace across the rounds). Also, it is common to set the mean of slopes to be 0 to avoid favoring any directions. Therefore we will set the mean of the prior for slope to 0.

priors_rslope_rate = c(
  prior(normal(1.82, 0.5), class = Intercept),
  prior(normal(0, 0.5), class = b),
  prior(normal(0, 0.5), class = sd),
  prior(lkj(2), class = cor))

priors_rslope_rate_zinb = c(
  prior(normal(1.82, 0.5), class = Intercept),
  prior(normal(0, 0.5), class = b),
  prior(normal(0, 0.5), class = sd),
  prior(lkj(2), class = cor),
  prior(normal(0, 0.5), class = zi), # on the logit scale
  prior(normal(0, 50), class = shape))


5.2.2 Model 3: condition x round (pair)

In the previous section, we analyzed the number of iconic gestures produced per trial. However, it is common to analyze the rate of iconic gestures per 100 words to account for the differences in the length of the trials and speech rate. Here, we will include the log of speech rate (number of words / 100) as an exposure variable and analyze the rate of iconic gestures per 100 words by condition. Note that the syntax for the exposure variable is different from the Poisson regression model; for negative binomial regression, the exposure variable is included with the rate() function.


5.2.2.1 Model Effect: round

nb_iconic_rate_cond_round = brm(num_iconic_gestures | rate(n_words_100) ~ 
                                  1 + condition * round + 
                                  (1+round|pair) + (1|target),
                                family = negbinomial(),
                                prior = priors_rslope_rate,
                                data = df_trial_info,
                                sample_prior = T,
                                warmup = nwu, iter = niter,
                                control = list(adapt_delta = 0.9, 
                                               max_treedepth = 15),
                                file = "models/nb_iconic_rate_cond_round")

nb_iconic_rate_cond_round_c = brm(num_iconic_gestures | rate(n_words_100) ~ 
                                    1 + condition * round_c + 
                                    (1+round_c|pair) + (1|target),
                                  family = negbinomial(),
                                  prior = priors_rslope_rate,
                                  data = df_trial_info,
                                  sample_prior = T,
                                  warmup = nwu, iter = niter,
                                  control = list(adapt_delta = 0.9, 
                                                 max_treedepth = 15),
                                  file = "models/nb_iconic_rate_cond_round_c")

nb_iconic_rate_cond_round_log = brm(num_iconic_gestures | rate(n_words_100) ~ 
                                      1 + condition * log_round_c + 
                                      (1+log_round_c|pair) + (1|target),
                                    family = negbinomial(),
                                    prior = priors_rslope_rate,
                                    data = df_trial_info,
                                    sample_prior = T,
                                    warmup = nwu, iter = niter,
                                    control = list(adapt_delta = 0.9, 
                                                   max_treedepth = 15),
                                    file = "models/nb_iconic_rate_cond_round_log")



### loo compare
if (!file.exists("models/loo_nb_iconic_rate_cond_round.rds")){
  nb_cond_round_loo = loo(nb_iconic_rate_cond_round)
  saveRDS(nb_cond_round_loo, file = "models/loo_nb_iconic_rate_cond_round.rds")
}

if (!file.exists("models/loo_nb_iconic_rate_cond_round_c.rds")){
  nb_cond_round_c_loo = loo(nb_iconic_rate_cond_round_c)
  saveRDS(nb_cond_round_c_loo, file = "models/loo_nb_iconic_rate_cond_round_c.rds")
}

if (!file.exists("models/loo_nb_iconic_rate_cond_round_log.rds")){
  nb_cond_round_log_loo = loo(nb_iconic_rate_cond_round_log)
  saveRDS(nb_cond_round_log_loo, file = "models/loo_nb_iconic_rate_cond_round_log.rds")
}

nb_cond_round_loo = readRDS("models/loo_nb_iconic_rate_cond_round.rds")
nb_cond_round_c_loo = readRDS("models/loo_nb_iconic_rate_cond_round_c.rds")
nb_cond_round_log_loo = readRDS("models/loo_nb_iconic_rate_cond_round_log.rds")

loo_compare(nb_cond_round_loo, nb_cond_round_c_loo, nb_cond_round_log_loo)
##                               elpd_diff se_diff
## nb_iconic_rate_cond_round_c     0.0       0.0  
## nb_iconic_rate_cond_round      -9.8       5.3  
## nb_iconic_rate_cond_round_log -13.0       4.2

The leave-one-out (LOO) Effect shows that centered round provides a larger predictive power (elpd_diff) and smaller standard error (se_diff) compared to the backward-difference coded round or centered log-transformed round. Therefore, we will use the centered round.


5.2.2.2 Model Effect: ZI or not

zinb_iconic_rate_cond_round_c = brm(num_iconic_gestures ~ 1 + condition * round_c + 
                                      offset(log(n_words_100)) +
                                      (1+round_c|pair) + (1|target),
                                    family = zero_inflated_negbinomial(),
                                    prior = priors_rslope_rate_zinb,
                                    data = df_trial_info,
                                    sample_prior = T,
                                    warmup = nwu, iter = niter,
                                    control = list(adapt_delta = 0.9, 
                                                   max_treedepth = 15),
                                    file = "models/zinb_iconic_rate_cond_round_c")



### loo compare
if (!file.exists("models/loo_zinb_iconic_rate_cond_round_c.rds")){
  zinb_cond_round_c_loo = loo(zinb_iconic_rate_cond_round_c)
  saveRDS(zinb_cond_round_c_loo, file = "models/loo_zinb_iconic_rate_cond_round_c.rds")
}

nb_cond_round_c_loo = readRDS("models/loo_nb_iconic_rate_cond_round_c.rds")
zinb_cond_round_c_loo = readRDS("models/loo_zinb_iconic_rate_cond_round_c.rds")

loo_compare(nb_cond_round_c_loo, zinb_cond_round_c_loo)
##                               elpd_diff se_diff
## zinb_iconic_rate_cond_round_c   0.0       0.0  
## nb_iconic_rate_cond_round_c   -13.7       6.2

The leave-one-out (LOO) Effect shows that zero-inflation model has a higher predictive power. As such, we will use zero-inflated negative binomial regression model.


5.2.2.3 Prior predictive check

zinb_iconic_rate_cond_round_c_prior = brm(num_iconic_gestures ~ 1 + condition * round_c + 
                                            offset(log(n_words_100)) +
                                            (1+round_c|pair) + (1|target),
                                          family = zero_inflated_negbinomial(),
                                          prior = priors_rslope_rate_zinb,
                                          sample_prior = "only",
                                          data = df_trial_info,
                                          file = "models/zinb_iconic_rate_cond_round_c_prior")

pp_check(zinb_iconic_rate_cond_round_c_prior, ndraws = 100, type = "bars") +
  coord_cartesian(xlim = c(0, 20),
                  ylim = c(0, 4000))

The prior predictive check shows that the model generates data that are somewhat similar to the observed data.


5.2.2.4 Fit the model

zinb_iconic_rate_cond_round_c = brm(num_iconic_gestures ~ 1 + condition * round_c + 
                                      offset(log(n_words_100)) +
                                      (1+round_c|pair) + (1|target),
                                    family = zero_inflated_negbinomial(),
                                    prior = priors_rslope_rate_zinb,
                                    data = df_trial_info,
                                    sample_prior = T,
                                    save_pars = save_pars(all = TRUE),
                                    warmup = nwu, iter = niter,
                                    control = list(adapt_delta = 0.9, 
                                                   max_treedepth = 15),
                                    file = "models/zinb_iconic_rate_cond_round_c")

model = zinb_iconic_rate_cond_round_c
summary(model)
##  Family: zero_inflated_negbinomial 
##   Links: mu = log; shape = identity; zi = identity 
## Formula: num_iconic_gestures ~ 1 + condition * round_c + offset(log(n_words_100)) + (1 + round_c | pair) + (1 | target) 
##    Data: df_trial_info (Number of observations: 4316) 
##   Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
##          total post-warmup draws = 72000
## 
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45) 
##                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)              0.66      0.07     0.53     0.82 1.00    18554
## sd(round_c)                0.11      0.02     0.08     0.15 1.00    26924
## cor(Intercept,round_c)     0.71      0.11     0.46     0.87 1.00    38762
##                        Tail_ESS
## sd(Intercept)             34315
## sd(round_c)               44094
## cor(Intercept,round_c)    50332
## 
## ~target (Number of levels: 16) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.10      0.03     0.06     0.16 1.00    26905    44443
## 
## Regression Coefficients:
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                     1.33      0.10     1.13     1.54 1.00    12630
## conditionAO_Asym              0.28      0.22    -0.15     0.70 1.00    14119
## conditionAsym_Sym            -0.06      0.22    -0.49     0.36 1.00    13682
## round_c                      -0.13      0.02    -0.17    -0.09 1.00    21391
## conditionAO_Asym:round_c      0.02      0.05    -0.07     0.11 1.00    24222
## conditionAsym_Sym:round_c    -0.03      0.04    -0.12     0.06 1.00    22801
##                           Tail_ESS
## Intercept                    23347
## conditionAO_Asym             25910
## conditionAsym_Sym            24753
## round_c                      36599
## conditionAO_Asym:round_c     40686
## conditionAsym_Sym:round_c    36775
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape    41.96     16.08    21.47    83.70 1.00   100817    52037
## zi        0.04      0.01     0.02     0.06 1.00    99853    44866
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayestestR::hdi(model)
# bayestestR::hdi(model, ci = 0.89)

The coefficients show that the condition does not have a significant effect on the rate of iconic gestures per 100 words. However, there was a significant decrease in the rate of iconic gestures per 100 words across the rounds. This means that the number of iconic gestures decreased more than the number of words did across the rounds. A formal hypothesis testing will be performed later using Bayes factor.


5.2.2.5 Visualize the posterior distributions

# p = plot_posterior(model, interaction = T)
# p


5.2.2.6 Hypothesis testing: Bayes factor

### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.1, 0.3, 1, 1.5)
bfs_cond_ao_asym = c()
bfs_cond_asym_sym = c()
bfs_round = c()

for (i in 1:length(prior_sd)){
  priors = c(
    prior(normal(1.82, 0.5), class = Intercept),
    set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b"),
    prior(normal(0, 0.5), class = sd),
    prior(lkj(2), class = cor),
    prior(normal(0, 0.5), class = zi), # on the logit scale
    prior(normal(0, 50), class = shape))
  
  fname = paste0("models/zinb_iconic_rate_cond_round_c_", prior_size[i])
  
  fit = brm(num_iconic_gestures ~ 
              1 + condition * round_c + 
              offset(log(n_words_100)) +
              (1+round_c|pair) + (1|target),
            family = zero_inflated_negbinomial(),
            prior = priors,
            data = df_trial_info,
            sample_prior = T,
            save_pars = save_pars(all = TRUE),
            warmup = nwu, iter = niter,
            control = list(adapt_delta = 0.9, 
                           max_treedepth = 15),
            file = fname)
  
  # BF for sym - asym
  h = hypothesis(fit, "conditionAO_Asym = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_cond_ao_asym = c(bfs_cond_ao_asym, bf)
  
  # BF for sym - ao
  h = hypothesis(fit, "conditionAsym_Sym = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_cond_asym_sym = c(bfs_cond_asym_sym, bf)
  
  # BF for round
  h = hypothesis(fit, "round_c = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_round = c(bfs_round, bf)
}

### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.5, prior_sd[3:4])

# BF for sym - asym
h = hypothesis(model, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym[3:5] = c(bf, bfs_cond_ao_asym[3:4])

# BF for sym - ao
h = hypothesis(model, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym[3:5] = c(bf, bfs_cond_asym_sym[3:4])

# BF for round
h = hypothesis(model, "round_c = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round[3:5] = c(bf, bfs_round[3:4])


### make a df for BFs
df_bf = data.frame(size = prior_size,
                   sd = prior_sd,
                   ao_asym = bfs_cond_ao_asym,
                   asym_sym = bfs_cond_asym_sym,
                   round = bfs_round) %>% 
  mutate(prior = paste0("N(0, ", sd, ")")) %>% 
  pivot_longer(cols = c("ao_asym", "asym_sym", "round"),
               names_to = "Effect",
               values_to = "BF10") %>% 
  mutate(Effect = factor(Effect,
                         levels = c("ao_asym", "asym_sym", "round"),
                         labels = c("AO-AsymAV", "AsymAV-SymAV", "Round")),
         Predictor = ifelse(Effect == "round", "Round", "Visibility"),
         BF10_log10 = log10(BF10))

df_bf %>% arrange(Effect, sd)
#### Plot BFs ####
ggplot(df_bf, aes(x = prior, y = BF10, group = Effect)) +
  geom_hline(yintercept = 1, linetype="dashed") +
  geom_point(aes(color=Effect)) +
  geom_line(aes(color=Effect)) +
  facet_wrap(vars(Predictor), scales="free_y") +
  theme_bw(base_size = 12)+
  theme(legend.position = "top")+
  scale_y_log10("Bayes factor (BF10)",
                breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
                labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
  xlab("prior")


5.2.2.7 Posterior predictive check

model = zinb_iconic_rate_cond_round_c
pp_check_sym = pp_check_each_condition(model, df_trial_info, "SymAV")
pp_check_asym = pp_check_each_condition(model, df_trial_info, "AsymAV")
pp_check_ao = pp_check_each_condition(model, df_trial_info, "AO")

gridExtra::grid.arrange(pp_check_sym, pp_check_asym, pp_check_ao, ncol = 2)

pp_check_r1 = pp_check_each_round(model, df_trial_info, "R1")
pp_check_r2 = pp_check_each_round(model, df_trial_info, "R2")
pp_check_r3 = pp_check_each_round(model, df_trial_info, "R3")
pp_check_r4 = pp_check_each_round(model, df_trial_info, "R4")
pp_check_r5 = pp_check_each_round(model, df_trial_info, "R5")
pp_check_r6 = pp_check_each_round(model, df_trial_info, "R6")

gridExtra::grid.arrange(pp_check_r1, pp_check_r2, pp_check_r3, pp_check_r4, pp_check_r5, pp_check_r6, ncol = 3)


5.2.2.8 Pair-wise Effect

emmeans(model, pairwise ~ condition)$contrasts
##  contrast       estimate lower.HPD upper.HPD
##  SymAV - AsymAV  -0.0633    -0.497     0.355
##  SymAV - AO       0.2134    -0.242     0.678
##  AsymAV - AO      0.2760    -0.140     0.708
## 
## Point estimate displayed: median 
## Results are given on the log (not the response) scale. 
## HPD interval probability: 0.95
emmeans(model, pairwise ~ condition, level = 0.89)$contrasts
##  contrast       estimate lower.HPD upper.HPD
##  SymAV - AsymAV  -0.0633    -0.405     0.286
##  SymAV - AO       0.2134    -0.163     0.584
##  AsymAV - AO      0.2760    -0.072     0.616
## 
## Point estimate displayed: median 
## Results are given on the log (not the response) scale. 
## HPD interval probability: 0.89


5.2.2.9 Visualize the model estimates

plot(conditional_effects(model), ask = FALSE)



6 =====Causal model of gestural alignment=====

Some experts in the field of statistics and causal inference have advised that researchers should build a causal model and examine which factors should be included and excluded from regression models if their aim is to infer the causal effects (e.g., McElreath, 2020, Pearl, 2010, Cinelli, Forney, & Pearl, 2020). Following this advice, we will build a causal model to infer the direct effect of speaker visibility on gestural alignment.

We assume the following causal model:

### causal model for gestural alignment rate
dag_gest <- dagitty('dag {
  visibility -> gest_align
  visibility -> n_words
  visibility -> n_gestures
  n_words -> lex_align
  n_words -> n_gestures
  round -> n_words
  round -> lex_align
  round -> gest_align
  round -> n_gestures
  n_gestures -> gest_align
  lex_align -> gest_align
}')

plot(dag_gest)

### check the minimam adjustment set
print("Direct effect of visibility on gestural alignment rate")
## [1] "Direct effect of visibility on gestural alignment rate"
adjustmentSets(dag_gest, exposure = "visibility", outcome = "gest_align", effect = "direct")
## { lex_align, n_gestures, round }
## { n_gestures, n_words, round }
print("Direct effect of lexical alignment rate on gestural alignment rate")
## [1] "Direct effect of lexical alignment rate on gestural alignment rate"
adjustmentSets(dag_gest, exposure = "lex_align", outcome = "gest_align", effect = "direct")
## { n_gestures, round, visibility }
## { n_words, round }
print("Direct effect of round on gestural alignment rate")
## [1] "Direct effect of round on gestural alignment rate"
adjustmentSets(dag_gest, exposure = "round", outcome = "gest_align", effect = "direct")
## { lex_align, n_gestures, visibility }

The minimum adjustment set for estimating the direct causal effect of speaker visibility on gestural alignment rate is {lex_align, n_words, round}. Note that because dagitty didn’t find any adjustment set for the direct effect of visibility on lexical alignment rate when we expected bidirectional causation between lexical and gestural alignment, we assumed a unidirecitonal causation from lexical alignment to gestural alignment only in this model. This will be reversed in the causal model for lexical alignment rate, such that we assume a unidirectional causation from gestural alignment to lexical alignment.

In addition, we are also interested in the direct effect of lexical alignment rate on gestural alignment rate. The minimum adjustment set for is {visibility, n_gestures, round}.

As the minimum adjustment sets for the direct effects of visibility, lexical alignment rate, and round on gestural alignment rate are identical (i.e., {visibility, lex_align, n_gestures, round}), we can estimate the direct effect of these variables on gestural alignment rate with the following model:

\[ gest\_align \: | \: \text{rate}(n\_iconic\_gesture) \sim \\ \text{visibility} * \text{round} + \text{n_lexical_alignment} + \\ (1+\text{round} | \text{participant}) + (1 | \text{item}) \]

Just to be thorough, we will first analyze the raw frequency and rate of gestural alignment. If you just want to see the analysis of the proportion of gestural alignment, see [-– Proportion of gestural alignment (n_gest_alignment / n_iconic)—].



7 =====Number of gestural alignment=====

7.1 —DataViz: number of gestural alignment—

7.1.1 bp: mean by condition

bp_mean_gest_alignment_by_cond = ggplot(data=trial_info_pair, 
                                        aes(x=condition, y=mean_num_gest_align, fill=condition)) +
  geom_jitter(aes(color = pair), 
              size = 1, alpha = 1, 
              width = 0.07, height = 0) +
  geom_boxplot(width = .4,
               outlier.shape = NA, alpha = 0.7) +
  geom_point(data = trial_info_cond, 
             aes(y = mean_num_gest_align), 
             size = 2, shape = 4) +
  scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
  labs(x = "Visibility",
       y = "Mean gestural alignment count") +
  theme_clean(base_size = 15) +
  theme(axis.text.x = element_text(colour = "black", size = 14),
        axis.text.y = element_text(colour = "black", size = 14),
        axis.title = element_text(size = 15, face = 'bold'),
        axis.title.x = element_text(vjust = -2),
        axis.title.y = element_text(vjust = 2),
        legend.position = "none",
        strip.text = element_text(size = 15, face = 'bold'),
        plot.background = element_blank(),
        plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))

ggplotly(bp_mean_gest_alignment_by_cond)


7.2 —Negative binomial regression—

7.2.1 Prior specification

We will set priors for the intercept based on the expected number of gestural alignment reported in Akamine et al. (2024). In the paper, the authors analyzed data from 19 dyads involving in the same task as the current study but in co-present interaction. The total number of gestural alignment reported was 1086, which was collected from 19 dyads, each performing 96 trials. This leads to the expected number of gestural alignment is 1086 / (19 * 96) = 0.6 (log(0.6) = -0.5). As such, we will set the prior for the intercept to be normal with a mean of -0.5 and a standard deviation of 0.5.

For the fixed effects, we will set unbiased priors with a mean of 0. We set different SDs for each effect because the scale of each predictor is different. For N. lexical alignment and N. iconic gestures, we set a prior of Normal(0, 0.2). This means that if the mean N. lex align or N. iconic gesture increase by 1, we expect the mean N. gestural alignment to change by -0.6 to 0.9, with the most likely size of change to be 0. As for the visibility condition and the round, we set a prior of Normal(0. 0.5).

  • N. lex align: Normal(0, 0.2) <– range of expected effect size in the original scale is [-0.6, 0.9]
  • N. iconic gestures: Normal(0, 0.2) <– [-0.6, 0.9]
  • visibility condition: Normal(0, 0.5) <– [-1.27, 3.13]
  • Round: Normal(0, 0.5) <– [-1.27, 3.13]

We will also specify a prior for the shape parameter to prevent the model from returning divergent transitions. We will set the prior to be normal with a mean of 0 and a wide standard deviation of 50.

For the standard deviation of the random effects, we set the prior to be normal with mean 0 and standard deviation 0.5. For the correlation between the random effects, we set the prior to be LKJ(2).

### pair
priors_rint_gest_align = c(
  prior(normal(0.6, 0.5), class = Intercept),
  prior(normal(0, 0.5), class = b),
  prior(normal(0, 0.2), class = b, coef = "num_lex_align"),
  prior(normal(0, 0.2), class = b, coef = "num_iconic_gestures"),
  prior(normal(0, 50), class = shape),
  prior(normal(0, 0.5), class = sd))

priors_rslope_gest_align = c(
  prior(normal(0.6, 0.5), class = Intercept),
  prior(normal(0, 0.5), class = b),
  prior(normal(0, 0.2), class = b, coef = "num_lex_align"),
  prior(normal(0, 0.2), class = b, coef = "num_iconic_gestures"),
  prior(normal(0, 50), class = shape),
  prior(normal(0, 0.5), class = sd),
  prior(lkj(2), class = cor))

priors_rslope_gest_align_zinb = c(
  prior(normal(0.6, 0.5), class = Intercept),
  prior(normal(0, 0.5), class = b),
  prior(normal(0, 0.2), class = b, coef = "num_lex_align"),
  prior(normal(0, 0.2), class = b, coef = "num_iconic_gestures"),
  prior(normal(0, 0.5), class = sd),
  prior(lkj(2), class = cor),
  prior(normal(0, 0.5), class = zi), # on the logit scale
  prior(normal(0, 50), class = shape))


7.2.2 Model 4: [pair] effect of lex align on gest align

7.2.2.1 Prior predictive check

zinb_gest_align_prior = brm(num_gestural_align ~
                              1 + condition + round + num_lex_align + num_iconic_gestures +
                              (1+round|pair) + (1|target),
                            family = zero_inflated_negbinomial(),
                            prior = priors_rslope_gest_align_zinb,
                            data = df_trial_info,
                            sample_prior = "only",
                            control = list(adapt_delta = 0.9, 
                                           max_treedepth = 20),
                            file = "models/zinb_gest_align_prior")

pp_check(zinb_gest_align_prior, ndraws = 100, type = "bars") +
  coord_cartesian(xlim = c(0, 20),
                  ylim = c(0, 4000))

The prior predictive check shows that the model expects fewer amount of 2, 3, and 4. This suggests that the zero-inflation prior may be too large or the mean for the intercept prior is too low. We will check the prior-posterior update plot and posterior predictive check to see if the model generates data that are similar to the observed data. If not, we will consider modifying the priors.


7.2.2.2 Fit the model

zinb_align_cond_round = brm(num_gestural_align ~ 
                              1 + condition + round + num_lex_align + num_iconic_gestures +
                              (1+round|pair) + (1|target),
                            family = zero_inflated_negbinomial(),
                            prior = priors_rslope_gest_align_zinb,
                            data = df_trial_info,
                            sample_prior = T,
                            save_pars = save_pars(all = TRUE),
                            warmup = nwu, iter = niter,
                            control = list(adapt_delta = 0.9, 
                                           max_treedepth = 15),
                            file = "models/zinb_align_cond_round")

model = zinb_align_cond_round
summary(model)
##  Family: zero_inflated_negbinomial 
##   Links: mu = log; shape = identity; zi = identity 
## Formula: num_gestural_align ~ 1 + condition + round + num_lex_align + num_iconic_gestures + (1 + round | pair) + (1 | target) 
##    Data: df_trial_info (Number of observations: 4316) 
##   Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
##          total post-warmup draws = 72000
## 
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45) 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)               1.04      0.12     0.82     1.30 1.00    13258
## sd(roundR12)                0.62      0.15     0.33     0.93 1.00    24974
## sd(roundR23)                0.29      0.12     0.05     0.52 1.00    14676
## sd(roundR34)                0.19      0.12     0.01     0.46 1.00    12288
## sd(roundR45)                0.22      0.13     0.01     0.49 1.00    16214
## sd(roundR56)                0.16      0.12     0.01     0.45 1.00    28127
## cor(Intercept,roundR12)     0.33      0.21    -0.12     0.69 1.00    38507
## cor(Intercept,roundR23)     0.47      0.25    -0.13     0.85 1.00    36116
## cor(roundR12,roundR23)      0.06      0.27    -0.47     0.59 1.00    43615
## cor(Intercept,roundR34)     0.18      0.31    -0.48     0.72 1.00    44648
## cor(roundR12,roundR34)      0.10      0.31    -0.52     0.66 1.00    47546
## cor(roundR23,roundR34)      0.16      0.31    -0.48     0.72 1.00    41117
## cor(Intercept,roundR45)     0.15      0.31    -0.48     0.69 1.00    50329
## cor(roundR12,roundR45)      0.16      0.30    -0.46     0.70 1.00    45069
## cor(roundR23,roundR45)      0.20      0.31    -0.46     0.74 1.00    31625
## cor(roundR34,roundR45)      0.11      0.32    -0.54     0.69 1.00    35440
## cor(Intercept,roundR56)    -0.08      0.33    -0.67     0.57 1.00    68954
## cor(roundR12,roundR56)      0.04      0.32    -0.58     0.64 1.00    62152
## cor(roundR23,roundR56)      0.04      0.33    -0.60     0.65 1.00    60123
## cor(roundR34,roundR56)      0.05      0.33    -0.59     0.66 1.00    50904
## cor(roundR45,roundR56)      0.03      0.33    -0.60     0.65 1.00    53360
##                         Tail_ESS
## sd(Intercept)              27474
## sd(roundR12)               27105
## sd(roundR23)               11704
## sd(roundR34)               20523
## sd(roundR45)               19835
## sd(roundR56)               30531
## cor(Intercept,roundR12)    46914
## cor(Intercept,roundR23)    32737
## cor(roundR12,roundR23)     47771
## cor(Intercept,roundR34)    48873
## cor(roundR12,roundR34)     53571
## cor(roundR23,roundR34)     49423
## cor(Intercept,roundR45)    48204
## cor(roundR12,roundR45)     49810
## cor(roundR23,roundR45)     44304
## cor(roundR34,roundR45)     48265
## cor(Intercept,roundR56)    53547
## cor(roundR12,roundR56)     53118
## cor(roundR23,roundR56)     57059
## cor(roundR34,roundR56)     56039
## cor(roundR45,roundR56)     56078
## 
## ~target (Number of levels: 16) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.22      0.06     0.13     0.35 1.00    21997    37924
## 
## Regression Coefficients:
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept              -1.93      0.17    -2.26    -1.58 1.00     7722    15921
## conditionAO_Asym        0.32      0.29    -0.25     0.88 1.00    13178    24373
## conditionAsym_Sym       0.25      0.27    -0.29     0.79 1.00    11227    21534
## roundR12                1.54      0.16     1.23     1.85 1.00    31455    40797
## roundR23               -0.10      0.11    -0.32     0.10 1.00    30196    40085
## roundR34               -0.25      0.11    -0.47    -0.06 1.00    35515    38965
## roundR45               -0.31      0.12    -0.56    -0.08 1.00    39319    44558
## roundR56               -0.18      0.12    -0.43     0.05 1.00    58580    50167
## num_lex_align           0.05      0.02     0.00     0.09 1.00    68008    56486
## num_iconic_gestures     0.20      0.01     0.18     0.23 1.00    41732    49580
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     3.65      0.53     2.77     4.83 1.00    49170    51339
## zi        0.01      0.01     0.00     0.03 1.00    72219    39837
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayestestR::hdi(model)

The coefficients show that the number of lexical alignment is a significant predictor of the number of gestural alignment.


7.2.2.3 Hypothesis testing: Bayes factor

### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.05, 0.1, 0.3, 0.5)
bfs_lex_align = c()

for (i in 1:length(prior_sd)){
  priors = c(
    prior(normal(0.6, 0.5), class = Intercept),
    prior(normal(0, 0.5), class = b),
    set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b", coef = "num_lex_align"),
    prior(normal(0, 0.2), class = b, coef = "num_iconic_gestures"),
    prior(normal(0, 0.5), class = sd),
    prior(lkj(2), class = cor),
    prior(normal(0, 0.5), class = zi),
    prior(normal(0, 50), class = shape))
  
  fname = paste0("models/zinb_align_cond_round_", prior_size[i])
  
  fit = brm(num_gestural_align ~ 
              1 + condition + round + num_lex_align + num_iconic_gestures +
              (1+round|pair) + (1|target),
            family = zero_inflated_negbinomial(),
            prior = priors,
            data = df_trial_info,
            sample_prior = T,
            save_pars = save_pars(all = TRUE),
            warmup = nwu, iter = niter,
            control = list(adapt_delta = 0.9, 
                           max_treedepth = 15),
            file = fname)
  
  ### BF for N. lex alignment
  h = hypothesis(fit, "num_lex_align = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_lex_align = c(bfs_lex_align, bf)
}


### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.2, prior_sd[3:4])

### BF for N. lex alignment
h = hypothesis(model, "num_lex_align = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align[3:5] = c(bf, bfs_lex_align[3:5])


### make a df for BFs
df_bf_lex = data.frame(size = prior_size,
                       sd = prior_sd,
                       lex_align = bfs_lex_align) %>% 
  mutate(prior = paste0("N(0, ", sd, ")")) %>% 
  pivot_longer(cols = c("lex_align"),
               names_to = "Effect",
               values_to = "BF10") %>% 
  mutate(Effect = "N. lex align",
         Predictor = "N. lex align")
#### Plot BFs ####
ggplot(filter(df_bf_lex),
       aes(x = factor(sd), y = BF10, group = Effect)) +
  geom_hline(yintercept = 1, linetype="dashed") +
  geom_point(aes(color=Effect)) +
  geom_line(aes(color=Effect)) +
  facet_wrap(vars(Predictor)) +
  theme_clean(base_size = 15) +
  theme(axis.text.x = element_text(colour = "black", size = 14),
        axis.text.y = element_text(colour = "black", size = 14),
        axis.title = element_text(size = 15, face = 'bold'),
        axis.title.x = element_text(vjust = -2),
        axis.title.y = element_text(vjust = 2),
        legend.position = "top",
        strip.text = element_text(size = 15, face = 'bold'),
        plot.background = element_blank(),
        plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
  scale_y_log10("Bayes factor (BF10)",
                breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
                labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
  xlab("SD for the prior")



8 =====<Chosen> Proportion of gestural alignment (n_gest_alignment / n_iconic)=====

We can analyze the proportion of gestural alignment in two ways: (1) modeling the rate of gestural alignment per iconic gesture using a negative binomial regression model and (2) modeling the probability of gestural alignment using a zero-one-inflated beta regression model.

Key differences in the two models are that the negative binomial regression model assumes that the rate of gestural alignment is a count variable, while the zero-one-inflated beta regression model assumes that the proportion of gestural alignment is a continuous variable bounded between 0 and 1. Also, while negative binomial regression models the rate of events, the zero-one-inflated beta regression models the probability of events. In the case of the proportion of gestural alignment, two models should yield similar results, but it is important to note that the two models are conceptually different.


8.1 —DataViz: proportion of gestural alignment—

8.1.1 bp: mean by condition

bp_mean_gest_alignment_prop_by_cond = ggplot(data=trial_info_pair, 
                                             aes(x=condition, y=mean_gest_align_prop, fill=condition)) +
  geom_jitter(aes(color = pair), 
              size = 1, alpha = 1, 
              width = 0.07, height = 0) +
  geom_boxplot(width = .3,
               outlier.shape = NA, alpha = 0.7) +
  geom_point(data = trial_info_cond, 
             aes(y = mean_gest_align_prop), 
             shape = 21, size = 3, fill = "white") +
  scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
  labs(x = "Visibility",
       y = "Mean gest align rate") +
  theme_clean(base_size = 15) +
  theme(axis.text.x = element_text(colour = "black", size = 14),
        axis.text.y = element_text(colour = "black", size = 14),
        axis.title = element_text(size = 15, face = 'bold'),
        axis.title.x = element_text(vjust = -2),
        axis.title.y = element_text(vjust = 2),
        legend.position = "none",
        strip.text = element_text(size = 15, face = 'bold'),
        plot.background = element_blank(),
        plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))

ggplotly(bp_mean_gest_alignment_prop_by_cond)


8.1.2 bp: mean by condition and round

pd = position_dodge(width = .75)

bp_mean_gest_alignment_prop_by_round_cond = 
  ggplot(data=trial_info_pair_round, 
         aes(x=round_n, y=mean_gest_align_prop, fill = condition)) +
  geom_jitter(aes(color = pair), 
              size = 0.5, alpha = 0.7, 
              width = 0.07, height = 0) +
  geom_boxplot(width = .5,
               outlier.shape = NA, alpha = 0.7) +
  geom_point(data = trial_info_cond_round, 
             aes(y = mean_gest_align_prop),
             position = pd,
             shape = 21, size = 2, fill = "white") +
  scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "Round",
       y = "Mean gest align rate") +
  theme_clean(base_size = 15) +
  theme(axis.text.x = element_text(colour = "black", size = 14),
        axis.text.y = element_text(colour = "black", size = 14),
        axis.title = element_text(size = 15, face = 'bold'),
        axis.title.x = element_text(vjust = -2),
        axis.title.y = element_text(vjust = 2),
        legend.position = "none",
        strip.text = element_text(size = 15, face = 'bold'),
        plot.background = element_blank(),
        plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
  facet_grid(cols = vars(condition))

ggplotly(bp_mean_gest_alignment_prop_by_round_cond)


8.2 Prepare data

To model the proportion of gestural alignment, we need to remove trials where the number of iconic gestures is 0. This is because dividing any numbers by 0 results in an undefined value (NA) and prevents model from running.

df_gest_align_posreg_prop = df_trial_info %>%
  filter(num_iconic_gestures > 0)

print(paste0("Number of rows before removing trials with no iconic gestures: ", nrow(df_trial_info)))
## [1] "Number of rows before removing trials with no iconic gestures: 4316"
print(paste0("Number of rows before after trials with no iconic gestures: ", nrow(df_gest_align_posreg_prop)))
## [1] "Number of rows before after trials with no iconic gestures: 2199"
print(paste0("Number of removed trials: ", nrow(df_trial_info) - nrow(df_gest_align_posreg_prop)))
## [1] "Number of removed trials: 2117"

8.3 Prior specification

We will set priors based on Akamine et al. (2024). As mentioned in the previous analysis, they detected 4413 iconic gestures and 1086 instances of gestural alignment. Dividing the number of gestural alignment by the number of iconic gestures gives 0.25 (1086/4413) (-1.39 on the log scale). This means that we expect 1 gestural alignment per 4 iconic gestures.

For the slope parameters, we set the mean of 0 with a SD of 0.5. This means that for example, we expect the mean difference between the AO and AsymAV conditions to range from -0.16 to 0.43.

### pair
priors_rint_gest_align_prop = c(
  prior(normal(-1.39, 0.5), class = Intercept),
  prior(normal(0, 0.5), class = b),
  prior(normal(0, 50), class = shape),
  prior(normal(0, 0.5), class = sd))

priors_rslope_gest_align_prop = c(
  prior(normal(-1.39, 0.5), class = Intercept),
  prior(normal(0, 0.5), class = b),
  prior(normal(0, 50), class = shape),
  prior(normal(0, 0.5), class = sd),
  prior(lkj(2), class = cor))

priors_rslope_gest_align_prop_zinb = c(
  prior(normal(-1.39, 0.5), class = Intercept),
  prior(normal(0, 0.5), class = b),
  prior(normal(0, 0.5), class = sd),
  prior(lkj(2), class = cor),
  prior(normal(0, 0.5), class = zi), # on the logit scale
  prior(normal(0, 50), class = shape))



8.4 —Negative binomial regression models—

8.4.1 Model 5: [pair] condition * round + num_lex_align

8.4.1.1 Model Effect: round

nb_align_rate_cond_round = brm(num_gestural_align | rate(num_iconic_gestures) ~ 
                                 1 + condition * round + num_lex_align +
                                 (1+round|pair) + (1|target),
                               family = negbinomial(),
                               prior = priors_rslope_gest_align_prop,
                               data = df_gest_align_posreg_prop,
                               sample_prior = T,
                               save_pars = save_pars(all = TRUE),
                               warmup = nwu, iter = niter,
                               control = list(adapt_delta = 0.9, 
                                              max_treedepth = 15),
                               file = "models/nb_align_rate_cond_round")

nb_align_rate_cond_round_c = brm(num_gestural_align | rate(num_iconic_gestures) ~ 
                                   1 + condition * round_c + num_lex_align +
                                   (1+round_c|pair) + (1|target),
                                 family = negbinomial(),
                                 prior = priors_rslope_gest_align_prop,
                                 data = df_gest_align_posreg_prop,
                                 sample_prior = T,
                                 save_pars = save_pars(all = TRUE),
                                 warmup = nwu, iter = niter,
                                 control = list(adapt_delta = 0.9, 
                                                max_treedepth = 15),
                                 file = "models/nb_align_rate_cond_round_c")

nb_align_rate_cond_round_log = brm(num_gestural_align | rate(num_iconic_gestures) ~ 
                                     1 + condition * log_round_c + num_lex_align +
                                     (1+log_round_c|pair) + (1|target),
                                   family = negbinomial(),
                                   prior = priors_rslope_gest_align_prop,
                                   data = df_gest_align_posreg_prop,
                                   sample_prior = T,
                                   save_pars = save_pars(all = TRUE),
                                   warmup = nwu, iter = niter,
                                   control = list(adapt_delta = 0.9, 
                                                  max_treedepth = 15),
                                   file = "models/nb_align_rate_cond_round_log")



### loo compare
if (!file.exists("models/loo_nb_align_rate_cond_round.rds")){
  nb_cond_round_loo = loo(nb_align_rate_cond_round)
  saveRDS(nb_cond_round_loo, file = "models/loo_nb_align_rate_cond_round.rds")
}

if (!file.exists("models/loo_nb_align_rate_cond_round_c.rds")){
  nb_cond_round_c_loo = loo(nb_align_rate_cond_round_c)
  saveRDS(nb_cond_round_c_loo, file = "models/loo_nb_align_rate_cond_round_c.rds")
}

if (!file.exists("models/loo_nb_align_rate_cond_round_log.rds")){
  nb_cond_round_log_loo = loo(nb_align_rate_cond_round_log)
  saveRDS(nb_cond_round_log_loo, file = "models/loo_nb_align_rate_cond_round_log.rds")
}

nb_cond_round_loo = readRDS("models/loo_nb_align_rate_cond_round.rds")
nb_cond_round_c_loo = readRDS("models/loo_nb_align_rate_cond_round_c.rds")
nb_cond_round_log_loo = readRDS("models/loo_nb_align_rate_cond_round_log.rds")

loo_compare(nb_cond_round_loo, nb_cond_round_c_loo, nb_cond_round_log_loo)
##                              elpd_diff se_diff
## nb_align_rate_cond_round        0.0       0.0 
## nb_align_rate_cond_round_log -131.4      16.3 
## nb_align_rate_cond_round_c   -238.7      22.7

The leave-one-out (LOO) Effect shows that the backward-difference coded round provides a far larger predictive power (elpd_diff) and a far smaller standard error (se_diff) compared to the centered round or centered log-transformed round. Thus, we will use the bd coded round as a predictor for further analyses.


8.4.1.2 Model Effect: ZI or not

zinb_align_rate_cond_round = brm(num_gestural_align ~ 
                                   1 + condition * round + num_lex_align +
                                   offset(log(num_iconic_gestures)) +
                                   (1+round|pair) + (1|target),
                                 family = zero_inflated_negbinomial(),
                                 prior = priors_rslope_gest_align_prop_zinb,
                                 data = df_gest_align_posreg_prop,
                                 sample_prior = T,
                                 warmup = nwu, iter = niter,
                                 control = list(adapt_delta = 0.9, 
                                                max_treedepth = 15),
                                 file = "models/zinb_align_rate_cond_round")



### loo compare
if (!file.exists("models/loo_zinb_align_rate_cond_round.rds")){
  zinb_cond_round_c_loo = loo(zinb_align_rate_cond_round)
  saveRDS(zinb_cond_round_c_loo, file = "models/loo_zinb_align_rate_cond_round.rds")
}

nb_cond_round_c_loo = readRDS("models/loo_nb_align_rate_cond_round.rds")
zinb_cond_round_c_loo = readRDS("models/loo_zinb_align_rate_cond_round.rds")

loo_compare(nb_cond_round_c_loo, zinb_cond_round_c_loo)
##                            elpd_diff se_diff
## nb_align_rate_cond_round    0.0       0.0   
## zinb_align_rate_cond_round -1.5       0.6

The loo comparsion shows that non-inflation model has a higher predictive power than the zero-inflation model. As such, we will use NB regression models for further analyses.


8.4.1.3 Prior predictive check

nb_gest_align_prop_prior = brm(num_gestural_align | rate(num_iconic_gestures) ~
                                 1 + condition + round + num_lex_align +
                                 (1+round|pair) + (1|target),
                               family = negbinomial(),
                               prior = priors_rslope_gest_align_prop,
                               data = df_gest_align_posreg_prop,
                               sample_prior = "only",
                               control = list(adapt_delta = 0.9, 
                                              max_treedepth = 20),
                               file = "models/nb_gest_align_prop_prior")

pp_check(nb_gest_align_prop_prior, ndraws = 100, type = "bars") +
  coord_cartesian(xlim = c(0, 20),
                  ylim = c(0, 3000))

The prior predictive check shows that the model expects fewer amount of 2, 3, and 4. This suggests that the zero-inflation prior may be too large or the mean for the intercept prior is too low. We will check the prior-posterior update plot and posterior predictive check to see if the model generates data that are similar to the observed data. If not, we will consider modifying the priors.


8.4.1.4 Fit the model

nb_align_rate_cond_round = brm(num_gestural_align | rate(num_iconic_gestures) ~ 
                                 1 + condition * round + num_lex_align +
                                 (1+round|pair) + (1|target),
                               family = negbinomial(),
                               prior = priors_rslope_gest_align_prop,
                               data = df_gest_align_posreg_prop,
                               sample_prior = T,
                               save_pars = save_pars(all = TRUE),
                               warmup = nwu, iter = niter,
                               control = list(adapt_delta = 0.9, 
                                              max_treedepth = 15),
                               file = "models/nb_align_rate_cond_round")

model = nb_align_rate_cond_round
summary(model)
##  Family: negbinomial 
##   Links: mu = log; shape = identity 
## Formula: num_gestural_align | rate(num_iconic_gestures) ~ 1 + condition * round + num_lex_align + (1 + round | pair) + (1 | target) 
##    Data: df_gest_align_posreg_prop (Number of observations: 2199) 
##   Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
##          total post-warmup draws = 72000
## 
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45) 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)               0.67      0.09     0.51     0.87 1.00    20129
## sd(roundR12)                0.64      0.12     0.43     0.90 1.00    32271
## sd(roundR23)                0.13      0.08     0.01     0.30 1.00    20697
## sd(roundR34)                0.08      0.07     0.00     0.25 1.00    31671
## sd(roundR45)                0.13      0.09     0.01     0.34 1.00    25468
## sd(roundR56)                0.13      0.10     0.00     0.36 1.00    35078
## cor(Intercept,roundR12)     0.03      0.21    -0.38     0.43 1.00    42101
## cor(Intercept,roundR23)     0.18      0.32    -0.49     0.74 1.00    63141
## cor(roundR12,roundR23)     -0.22      0.30    -0.73     0.44 1.00    56062
## cor(Intercept,roundR34)    -0.00      0.33    -0.63     0.63 1.00    86753
## cor(roundR12,roundR34)     -0.09      0.32    -0.67     0.56 1.00    75896
## cor(roundR23,roundR34)      0.01      0.33    -0.63     0.63 1.00    65940
## cor(Intercept,roundR45)     0.04      0.33    -0.60     0.65 1.00    79946
## cor(roundR12,roundR45)      0.09      0.31    -0.53     0.65 1.00    71978
## cor(roundR23,roundR45)      0.01      0.33    -0.62     0.64 1.00    57760
## cor(roundR34,roundR45)     -0.04      0.34    -0.66     0.61 1.00    53393
## cor(Intercept,roundR56)    -0.06      0.34    -0.68     0.60 1.00    85023
## cor(roundR12,roundR56)      0.05      0.32    -0.57     0.64 1.00    79246
## cor(roundR23,roundR56)      0.00      0.33    -0.63     0.63 1.00    65196
## cor(roundR34,roundR56)     -0.01      0.33    -0.64     0.63 1.00    58872
## cor(roundR45,roundR56)     -0.02      0.33    -0.65     0.62 1.00    56595
##                         Tail_ESS
## sd(Intercept)              38159
## sd(roundR12)               47486
## sd(roundR23)               26398
## sd(roundR34)               34620
## sd(roundR45)               31341
## sd(roundR56)               32642
## cor(Intercept,roundR12)    50684
## cor(Intercept,roundR23)    50618
## cor(roundR12,roundR23)     49539
## cor(Intercept,roundR34)    52609
## cor(roundR12,roundR34)     53428
## cor(roundR23,roundR34)     57498
## cor(Intercept,roundR45)    53014
## cor(roundR12,roundR45)     54864
## cor(roundR23,roundR45)     57810
## cor(roundR34,roundR45)     56422
## cor(Intercept,roundR56)    53281
## cor(roundR12,roundR56)     54104
## cor(roundR23,roundR56)     56138
## cor(roundR34,roundR56)     58815
## cor(roundR45,roundR56)     58976
## 
## ~target (Number of levels: 16) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.06      0.04     0.00     0.14 1.00    22309    20923
## 
## Regression Coefficients:
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                     -1.55      0.11    -1.77    -1.33 1.00    13026
## conditionAO_Asym               0.38      0.23    -0.08     0.84 1.00    16239
## conditionAsym_Sym              0.39      0.22    -0.06     0.83 1.00    12812
## roundR12                       1.61      0.14     1.33     1.90 1.00    38765
## roundR23                       0.24      0.08     0.09     0.39 1.00    54571
## roundR34                       0.02      0.08    -0.14     0.19 1.00    58090
## roundR45                      -0.12      0.10    -0.32     0.08 1.00    54198
## roundR56                      -0.04      0.12    -0.27     0.19 1.00    64843
## num_lex_align                 -0.01      0.01    -0.04     0.02 1.00   100913
## conditionAO_Asym:roundR12     -0.07      0.28    -0.62     0.49 1.00    43120
## conditionAsym_Sym:roundR12     0.13      0.26    -0.39     0.64 1.00    36961
## conditionAO_Asym:roundR23      0.01      0.17    -0.32     0.34 1.00    57584
## conditionAsym_Sym:roundR23    -0.01      0.15    -0.30     0.29 1.00    55291
## conditionAO_Asym:roundR34      0.04      0.19    -0.33     0.40 1.00    55087
## conditionAsym_Sym:roundR34     0.15      0.16    -0.17     0.47 1.00    54834
## conditionAO_Asym:roundR45      0.24      0.22    -0.18     0.67 1.00    59488
## conditionAsym_Sym:roundR45    -0.14      0.19    -0.51     0.23 1.00    58188
## conditionAO_Asym:roundR56      0.19      0.26    -0.32     0.70 1.00    70505
## conditionAsym_Sym:roundR56     0.11      0.21    -0.30     0.51 1.00    72689
##                            Tail_ESS
## Intercept                     26179
## conditionAO_Asym              30017
## conditionAsym_Sym             24995
## roundR12                      48607
## roundR23                      50032
## roundR34                      50833
## roundR45                      51064
## roundR56                      54866
## num_lex_align                 55847
## conditionAO_Asym:roundR12     52077
## conditionAsym_Sym:roundR12    47901
## conditionAO_Asym:roundR23     55257
## conditionAsym_Sym:roundR23    51500
## conditionAO_Asym:roundR34     55384
## conditionAsym_Sym:roundR34    54773
## conditionAO_Asym:roundR45     53403
## conditionAsym_Sym:roundR45    54752
## conditionAO_Asym:roundR56     55218
## conditionAsym_Sym:roundR56    57241
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape    77.70     29.04    31.65   143.21 1.00    96781    51091
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayestestR::hdi(model)
# bayestestR::hdi(model, ci = 0.89)

The coefficients show that the proportion of gestural alignment was significantly higher in the SymAV condition than in the AO condition. Also, the rate of gestural alignment significantly increased from R1–R2 and R2–R3 and stabilized afterwards.


8.4.1.5 Hypothesis testing: Bayes factor

### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.1, 0.3, 0.7, 1)
bfs_cond_ao_asym = c()
bfs_cond_asym_sym = c()
bfs_round12 = c()
bfs_round23 = c()
bfs_round34 = c()
bfs_round45 = c()
bfs_round56 = c()
bfs_lex_align = c()

for (i in 1:length(prior_sd)){
  priors = c(
    prior(normal(-1.39, 0.5), class = Intercept),
    set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b"),
    prior(normal(0, 0.5), class = sd),
    prior(lkj(2), class = cor),
    prior(normal(0, 50), class = shape))
  
  fname = paste0("models/nb_align_rate_cond_round_", prior_size[i])
  
  fit = brm(num_gestural_align | rate(num_iconic_gestures) ~ 
              1 + condition * round + num_lex_align +
              (1+round|pair) + (1|target),
            family = negbinomial(),
            prior = priors,
            data = df_gest_align_posreg_prop,
            sample_prior = T,
            save_pars = save_pars(all = TRUE),
            warmup = nwu, iter = niter,
            control = list(adapt_delta = 0.9, 
                           max_treedepth = 15),
            file = fname)
  
  ### BF for visibility conditions
  # sym - asym
  h = hypothesis(fit, "conditionAO_Asym = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_cond_ao_asym = c(bfs_cond_ao_asym, bf)
  
  # sym - ao
  h = hypothesis(fit, "conditionAsym_Sym = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_cond_asym_sym = c(bfs_cond_asym_sym, bf)
  
  ### BF for rounds
  # R1 - R2
  h = hypothesis(model, "roundR12 = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_round12 = c(bfs_round12, bf)
  
  # R2 - R3
  h = hypothesis(model, "roundR23 = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_round23 = c(bfs_round23, bf)
  
  # R3 - R4
  h = hypothesis(model, "roundR34 = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_round34 = c(bfs_round34, bf)
  
  # R4 - R5
  h = hypothesis(model, "roundR45 = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_round45 = c(bfs_round45, bf)
  
  # R5 - R6
  h = hypothesis(model, "roundR56 = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_round56 = c(bfs_round56, bf)
  
  ### BF for N. lex alignment
  h = hypothesis(fit, "num_lex_align = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_lex_align = c(bfs_lex_align, bf)
}


### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.5, prior_sd[3:4])

### BF for visibility
# sym - asym
h = hypothesis(model, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym[3:5] = c(bf, bfs_cond_ao_asym[3:4])

# sym - ao
h = hypothesis(model, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym[3:5] = c(bf, bfs_cond_asym_sym[3:4])

### BF for rounds
# R1 - R2
h = hypothesis(model, "roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round12[3:5] = c(bf, bfs_round12[3:4])

# R2 - R3
h = hypothesis(model, "roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round23[3:5] = c(bf, bfs_round23[3:4])

# R3 - R4
h = hypothesis(model, "roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round34[3:5] = c(bf, bfs_round34[3:4])

# R4 - R5
h = hypothesis(model, "roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round45[3:5] = c(bf, bfs_round45[3:4])

# R5 - R6
h = hypothesis(model, "roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round56[3:5] = c(bf, bfs_round56[3:4])

### BF for N. lex alignment
h = hypothesis(model, "num_lex_align = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align[3:5] = c(bf, bfs_lex_align[3:5])


### make a df for BFs
df_bf = data.frame(size = prior_size,
                   sd = prior_sd,
                   ao_asym = bfs_cond_ao_asym,
                   asym_sym = bfs_cond_asym_sym,
                   round12 = bfs_round12,
                   round23 = bfs_round23,
                   round34 = bfs_round34,
                   round45 = bfs_round45,
                   round56 = bfs_round56,
                   lex_align = bfs_lex_align) %>% 
  mutate(prior = paste0("N(0, ", sd, ")")) %>% 
  pivot_longer(cols = c("ao_asym", "asym_sym", 
                        "round12", "round23", "round34", "round45", "round56",
                        "lex_align"),
               names_to = "Effect",
               values_to = "BF10") %>% 
  mutate(Predictor = ifelse(grepl("round", Effect), "Round", 
                            ifelse(Effect == "lex_align", "N. lex align", 
                                   "Visibility")))

df_bf$Effect = recode(df_bf$Effect,
                      ao_asym = "AO--AsymAV",
                      asym_sym = "AsymAV--SymAV",
                      round12 = "R1--R2",
                      round23 = "R2--R3",
                      round34 = "R3--R4",
                      round45 = "R4--R5",
                      round56 = "R5--R6",
                      lex_align = "N. lex align")
#### Plot BFs ####
ggplot(filter(df_bf, Effect!="R1--R2"), #exclude R1--R2 because BF is too huge
       aes(x = factor(sd), y = BF10, group = Effect)) +
  geom_hline(yintercept = 1, linetype="dashed") +
  geom_point(aes(color=Effect)) +
  geom_line(aes(color=Effect)) +
  facet_wrap(vars(Predictor)) +
  theme_clean(base_size = 15) +
  theme(axis.text.x = element_text(colour = "black", size = 14),
        axis.text.y = element_text(colour = "black", size = 14),
        axis.title = element_text(size = 15, face = 'bold'),
        axis.title.x = element_text(vjust = -2),
        axis.title.y = element_text(vjust = 2),
        legend.position = "top",
        strip.text = element_text(size = 15, face = 'bold'),
        plot.background = element_blank(),
        plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
  scale_y_log10("Bayes factor (BF10)",
                breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
                labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
  xlab("SD for the prior")


8.4.1.6 Probability of direction

p_direction(model)


8.4.1.7 Pair-wise effect

emmeans(model, pairwise ~ condition)$contrasts
##  contrast       estimate lower.HPD upper.HPD
##  SymAV - AsymAV    0.388   -0.0568     0.829
##  SymAV - AO        0.770    0.2787     1.265
##  AsymAV - AO       0.382   -0.0787     0.832
## 
## Results are averaged over the levels of: round 
## Point estimate displayed: median 
## Results are given on the log (not the response) scale. 
## HPD interval probability: 0.95
emmeans(model, pairwise ~ condition, level = 0.89)$contrasts
##  contrast       estimate lower.HPD upper.HPD
##  SymAV - AsymAV    0.388     0.028     0.744
##  SymAV - AO        0.770     0.382     1.182
##  AsymAV - AO       0.382     0.012     0.750
## 
## Results are averaged over the levels of: round 
## Point estimate displayed: median 
## Results are given on the log (not the response) scale. 
## HPD interval probability: 0.89


8.4.1.8 Prior-posterior update plot

post_sample = as_draws_df(model)
pp_update_plot(post_sample, model_type="nb")


8.4.1.9 Posterior predictive check

pp_check_overall = pp_check(model, ndraws = 100, type = "bars") +
  coord_cartesian(xlim = c(0, 20),
                  ylim = c(0, 3000))
pp_check_sym = pp_check_each_condition(model, df_gest_align_posreg_prop, "SymAV")
pp_check_asym = pp_check_each_condition(model, df_gest_align_posreg_prop, "AsymAV")
pp_check_ao = pp_check_each_condition(model, df_gest_align_posreg_prop, "AO")

gridExtra::grid.arrange(pp_check_overall, pp_check_sym, 
                        pp_check_asym, pp_check_ao, 
                        ncol = 2)

Although the model prediction is not perfect, this model had a higher predictive power than the negative binomial model. As such, we will use this model.


8.4.1.10 Visualize the model estimates

plot(conditional_effects(model), ask=FALSE)



8.4.2 Model 6: [pair] condition_sum * round + num_lex_align

Here, we will run another model to compare the rate of gestural alignment between the SymAV and AO conditions.


8.4.2.1 Prior predictive check

nb_gest_align_prop_sum_prior = brm(num_gestural_align | rate(num_iconic_gestures) ~
                                     1 + condition_sum + round + num_lex_align +
                                     (1+round|pair) + (1|target),
                                   family = negbinomial(),
                                   prior = priors_rslope_gest_align_prop,
                                   data = df_gest_align_posreg_prop,
                                   sample_prior = "only",
                                   control = list(adapt_delta = 0.9, 
                                                  max_treedepth = 20),
                                   file = "models/nb_gest_align_prop_sum_prior")

pp_check(nb_gest_align_prop_sum_prior, ndraws = 100, type = "bars") +
  coord_cartesian(xlim = c(0, 20),
                  ylim = c(0, 3000))

The prior predictive check shows that the model expects fewer amount of 2, 3, and 4. This suggests that the zero-inflation prior may be too large or the mean for the intercept prior is too low. We will check the prior-posterior update plot and posterior predictive check to see if the model generates data that are similar to the observed data. If not, we will consider modifying the priors.


8.4.2.2 Fit the model

nb_align_rate_cond_round_sum = brm(num_gestural_align | rate(num_iconic_gestures) ~ 
                                     1 + condition_sum * round + num_lex_align +
                                     (1+round|pair) + (1|target),
                                   family = negbinomial(),
                                   prior = priors_rslope_gest_align_prop,
                                   data = df_gest_align_posreg_prop,
                                   sample_prior = T,
                                   save_pars = save_pars(all = TRUE),
                                   warmup = nwu, iter = niter,
                                   control = list(adapt_delta = 0.9, 
                                                  max_treedepth = 15),
                                   file = "models/nb_align_rate_cond_round_sum")

model = nb_align_rate_cond_round_sum
summary(model)
##  Family: negbinomial 
##   Links: mu = log; shape = identity 
## Formula: num_gestural_align | rate(num_iconic_gestures) ~ 1 + condition_sum * round + num_lex_align + (1 + round | pair) + (1 | target) 
##    Data: df_gest_align_posreg_prop (Number of observations: 2199) 
##   Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
##          total post-warmup draws = 72000
## 
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45) 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)               0.67      0.09     0.51     0.87 1.00    18575
## sd(roundR12)                0.64      0.12     0.43     0.90 1.00    31125
## sd(roundR23)                0.12      0.08     0.01     0.30 1.00    20919
## sd(roundR34)                0.08      0.07     0.00     0.24 1.00    32022
## sd(roundR45)                0.12      0.09     0.00     0.33 1.00    26026
## sd(roundR56)                0.13      0.10     0.00     0.37 1.00    34132
## cor(Intercept,roundR12)     0.04      0.21    -0.38     0.44 1.00    40280
## cor(Intercept,roundR23)     0.18      0.32    -0.49     0.73 1.00    65536
## cor(roundR12,roundR23)     -0.22      0.31    -0.73     0.44 1.00    59711
## cor(Intercept,roundR34)    -0.01      0.34    -0.64     0.63 1.00    93535
## cor(roundR12,roundR34)     -0.09      0.32    -0.67     0.56 1.00    78483
## cor(roundR23,roundR34)      0.01      0.33    -0.62     0.64 1.00    67466
## cor(Intercept,roundR45)     0.04      0.33    -0.60     0.65 1.00    80878
## cor(roundR12,roundR45)      0.09      0.31    -0.54     0.66 1.00    78177
## cor(roundR23,roundR45)      0.01      0.33    -0.62     0.64 1.00    57583
## cor(roundR34,roundR45)     -0.04      0.33    -0.66     0.61 1.00    54005
## cor(Intercept,roundR56)    -0.06      0.34    -0.68     0.59 1.00    93756
## cor(roundR12,roundR56)      0.05      0.32    -0.57     0.64 1.00    85093
## cor(roundR23,roundR56)     -0.00      0.33    -0.63     0.63 1.00    67602
## cor(roundR34,roundR56)     -0.01      0.33    -0.63     0.62 1.00    60265
## cor(roundR45,roundR56)     -0.02      0.33    -0.64     0.62 1.00    58058
##                         Tail_ESS
## sd(Intercept)              34700
## sd(roundR12)               48295
## sd(roundR23)               26610
## sd(roundR34)               31291
## sd(roundR45)               32316
## sd(roundR56)               33988
## cor(Intercept,roundR12)    50279
## cor(Intercept,roundR23)    51923
## cor(roundR12,roundR23)     48590
## cor(Intercept,roundR34)    54422
## cor(roundR12,roundR34)     54685
## cor(roundR23,roundR34)     57533
## cor(Intercept,roundR45)    51517
## cor(roundR12,roundR45)     55396
## cor(roundR23,roundR45)     57358
## cor(roundR34,roundR45)     57922
## cor(Intercept,roundR56)    54259
## cor(roundR12,roundR56)     55731
## cor(roundR23,roundR56)     56162
## cor(roundR34,roundR56)     59440
## cor(roundR45,roundR56)     59143
## 
## ~target (Number of levels: 16) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.06      0.03     0.00     0.14 1.00    21296    20498
## 
## Regression Coefficients:
##                                Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                         -1.54      0.11    -1.77    -1.33 1.00
## condition_sumAO_Sym                0.65      0.23     0.19     1.11 1.00
## condition_sumAsym_Sym              0.27      0.23    -0.18     0.72 1.00
## roundR12                           1.62      0.14     1.33     1.90 1.00
## roundR23                           0.24      0.08     0.09     0.39 1.00
## roundR34                           0.03      0.08    -0.13     0.19 1.00
## roundR45                          -0.12      0.10    -0.32     0.08 1.00
## roundR56                          -0.03      0.12    -0.27     0.20 1.00
## num_lex_align                     -0.01      0.01    -0.04     0.02 1.00
## condition_sumAO_Sym:roundR12       0.03      0.28    -0.52     0.57 1.00
## condition_sumAsym_Sym:roundR12     0.13      0.26    -0.39     0.65 1.00
## condition_sumAO_Sym:roundR23       0.01      0.17    -0.32     0.34 1.00
## condition_sumAsym_Sym:roundR23    -0.01      0.15    -0.30     0.29 1.00
## condition_sumAO_Sym:roundR34       0.16      0.18    -0.20     0.52 1.00
## condition_sumAsym_Sym:roundR34     0.15      0.16    -0.17     0.47 1.00
## condition_sumAO_Sym:roundR45       0.13      0.22    -0.29     0.56 1.00
## condition_sumAsym_Sym:roundR45    -0.15      0.19    -0.52     0.22 1.00
## condition_sumAO_Sym:roundR56       0.25      0.26    -0.26     0.76 1.00
## condition_sumAsym_Sym:roundR56     0.08      0.21    -0.33     0.49 1.00
##                                Bulk_ESS Tail_ESS
## Intercept                         13062    24757
## condition_sumAO_Sym               13128    25189
## condition_sumAsym_Sym             11201    23088
## roundR12                          38899    48092
## roundR23                          54168    50887
## roundR34                          61182    52654
## roundR45                          58074    53364
## roundR56                          66159    55148
## num_lex_align                     95708    56583
## condition_sumAO_Sym:roundR12      43962    53222
## condition_sumAsym_Sym:roundR12    36727    47437
## condition_sumAO_Sym:roundR23      58050    53368
## condition_sumAsym_Sym:roundR23    60204    55326
## condition_sumAO_Sym:roundR34      60163    56713
## condition_sumAsym_Sym:roundR34    60303    55921
## condition_sumAO_Sym:roundR45      60850    56172
## condition_sumAsym_Sym:roundR45    60868    55720
## condition_sumAO_Sym:roundR56      70417    57299
## condition_sumAsym_Sym:roundR56    74686    56413
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape    77.56     29.03    31.62   143.29 1.00   107322    49549
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayestestR::hdi(model)
# bayestestR::hdi(model, ci = 0.89)

The coefficients show that the proportion of gestural alignment was significantly higher in the SymAV condition than in the AO condition. Also, the rate of gestural alignment significantly increased from R1–R2 and R2–R3 and stabilized afterwards.


8.4.2.3 Visualize posterior distributions

main_model = nb_align_rate_cond_round
model1 = zinb_align_cond_round
plot_posterior(main_model, model1, model)


8.4.2.4 Hypothesis testing: Bayes factor

### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.1, 0.3, 0.7, 1)
bfs_cond_ao_sym = c()


for (i in 1:length(prior_sd)){
  priors = c(
    prior(normal(-1.39, 0.5), class = Intercept),
    set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b"),
    prior(normal(0, 0.5), class = sd),
    prior(lkj(2), class = cor),
    prior(normal(0, 50), class = shape))
  
  fname = paste0("models/nb_align_rate_cond_round_sum_", prior_size[i])
  
  fit = brm(num_gestural_align | rate(num_iconic_gestures) ~ 
              1 + condition_sum * round + num_lex_align +
              (1+round|pair) + (1|target),
            family = negbinomial(),
            prior = priors,
            data = df_gest_align_posreg_prop,
            sample_prior = T,
            save_pars = save_pars(all = TRUE),
            warmup = nwu, iter = niter,
            control = list(adapt_delta = 0.9, 
                           max_treedepth = 15),
            file = fname)
  
  ### BF for visibility conditions
  # BF for sym - ao
  h = hypothesis(fit, "condition_sumAO_Sym = 0")
  bf = 1 / abs(h$hypothesis$Evid.Ratio)
  bfs_cond_ao_sym = c(bfs_cond_ao_sym, bf)
}

### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.5, prior_sd[3:4])

### BF for visibility conditions
# sym - ao
h = hypothesis(model, "condition_sumAO_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_sym[3:5] = c(bf, bfs_cond_ao_sym[3:4])


### make a df for BFs
df_bf_temp = data.frame(size = prior_size,
                        sd = prior_sd,
                        Effect = "AO--SymAV",
                        BF10 = bfs_cond_ao_sym) %>% 
  mutate(prior = paste0("N(0, ", sd, ")"),
         Predictor = "Visibility")

df_bf_new = df_bf %>% 
  filter(Effect != "N. lex align") %>% 
  rbind(df_bf_temp) %>% 
  rbind(df_bf_lex) %>% 
  mutate(Effect = factor(Effect,
                         levels = c("AO--SymAV", "AO--AsymAV", "AsymAV--SymAV", 
                                    "R1--R2", "R2--R3", "R3--R4", "R4--R5", "R5--R6",
                                    "N. lex align")),
         Predictor = factor(Predictor,
                            levels = c("Visibility", "Round", "N. lex align")))

df_bf_new %>% arrange(Effect, sd)
#### Plot BFs ####
ggplot(filter(df_bf_new, Effect!="R1--R2"), #exclude R1--R2 because BF is too huge
       aes(x = factor(sd), y = BF10, group = Effect)) +
  geom_hline(yintercept = 1, linetype="dashed") +
  geom_point(aes(color=Effect)) +
  geom_line(aes(color=Effect)) +
  facet_wrap(vars(Predictor), scales = "free_x") +
  theme_clean(base_size = 15) +
  theme(axis.text.x = element_text(colour = "black", size = 14),
        axis.text.y = element_text(colour = "black", size = 14),
        axis.title = element_text(size = 15, face = 'bold'),
        axis.title.x = element_text(vjust = -2),
        axis.title.y = element_text(vjust = 2),
        legend.position = "top",
        legend.title=element_text(size=14), 
        legend.text=element_text(size=13),
        strip.text = element_text(size = 15, face = 'bold'),
        plot.background = element_blank(),
        plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
  scale_y_log10("Bayes factor (BF10)",
                breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
                labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
  xlab("SD for the prior")


8.4.2.5 Probability of direction

p_direction(model)


8.4.2.6 Prior-posterior update plot

post_sample = as_draws_df(model)
pp_update_plot(post_sample, model_type="nb")


8.4.2.7 Posterior predictive check

pp_check_overall = pp_check(model, ndraws = 100, type = "bars") +
  coord_cartesian(xlim = c(0, 20),
                  ylim = c(0, 3000))
pp_check_sym = pp_check_each_condition(model, df_gest_align_posreg_prop, "SymAV")
pp_check_asym = pp_check_each_condition(model, df_gest_align_posreg_prop, "AsymAV")
pp_check_ao = pp_check_each_condition(model, df_gest_align_posreg_prop, "AO")

gridExtra::grid.arrange(pp_check_overall, pp_check_sym, 
                        pp_check_asym, pp_check_ao, 
                        ncol = 2)

Although the model prediction is not perfect, this model had a higher predictive power than the negative binomial model. As such, we will use this model.


8.4.2.8 Visualize the model estimates

plot(conditional_effects(model), ask=FALSE)



9 =====Correlation btw lexical and gestural alignment=====

cor = pcor.test(df_trial_info$num_lex_align, df_trial_info$num_gestural_align, df_trial_info$log_round_c)
cor



10 =====Session info=====

sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22631)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] rstan_2.32.6       StanHeaders_2.32.7 ppcor_1.1          MASS_7.3-58.1     
##  [5] doParallel_1.0.17  iterators_1.0.14   foreach_1.5.2      emmeans_1.10.7    
##  [9] tidybayes_3.0.7    bayestestR_0.15.2  brms_2.22.0        Rcpp_1.0.12       
## [13] plotly_4.10.4      rsvg_2.6.2         DiagrammeRsvg_0.1  DiagrammeR_1.0.11 
## [17] dagitty_0.3-4      ggh4x_0.3.1        ggthemes_5.1.0     hypr_0.2.8        
## [21] plotrix_3.8-4      lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1     
## [25] dplyr_1.1.4        purrr_1.0.2        readr_2.1.5        tidyr_1.3.1       
## [29] tibble_3.2.1       ggplot2_3.5.2      tidyverse_2.0.0   
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.1-0     estimability_1.5.1   QuickJSR_1.1.3      
##  [4] rstudioapi_0.17.1    farver_2.1.1         svUnit_1.0.6        
##  [7] bit64_4.0.5          mvtnorm_1.2-4        bridgesampling_1.1-2
## [10] codetools_0.2-18     cachem_1.0.8         knitr_1.49          
## [13] bayesplot_1.11.1     jsonlite_1.8.8       ggdist_3.3.2        
## [16] compiler_4.2.2       httr_1.4.7           backports_1.4.1     
## [19] Matrix_1.5-1         fastmap_1.1.1        lazyeval_0.2.2      
## [22] cli_3.6.2            visNetwork_2.1.2     htmltools_0.5.8.1   
## [25] tools_4.2.2          coda_0.19-4.1        gtable_0.3.6        
## [28] glue_1.7.0           reshape2_1.4.4       posterior_1.6.1     
## [31] V8_6.0.6             jquerylib_0.1.4      vctrs_0.6.5         
## [34] nlme_3.1-160         crosstalk_1.2.1      insight_1.1.0       
## [37] tensorA_0.36.2.1     xfun_0.53            ps_1.7.6            
## [40] timechange_0.3.0     lifecycle_1.0.4      scales_1.3.0        
## [43] vroom_1.6.5          ragg_1.3.0           hms_1.1.3           
## [46] Brobdingnag_1.2-9    inline_0.3.21        RColorBrewer_1.1-3  
## [49] yaml_2.3.8           curl_6.2.1           gridExtra_2.3       
## [52] loo_2.8.0            sass_0.4.9           stringi_1.8.3       
## [55] checkmate_2.3.1      pkgbuild_1.4.6       boot_1.3-28         
## [58] cmdstanr_0.8.1       rlang_1.1.3          pkgconfig_2.0.3     
## [61] systemfonts_1.0.6    matrixStats_1.3.0    distributional_0.5.0
## [64] pracma_2.4.4         evaluate_1.0.3       lattice_0.20-45     
## [67] rstantools_2.4.0     htmlwidgets_1.6.4    labeling_0.4.3      
## [70] processx_3.8.4       bit_4.0.5            tidyselect_1.2.1    
## [73] plyr_1.8.9           magrittr_2.0.3       R6_2.6.1            
## [76] generics_0.1.3       pillar_1.10.1        withr_3.0.2         
## [79] datawizard_1.0.1     abind_1.4-8          crayon_1.5.3        
## [82] arrayhelpers_1.1-0   tzdb_0.4.0           rmarkdown_2.29      
## [85] grid_4.2.2           data.table_1.15.4    digest_0.6.35       
## [88] xtable_1.8-4         textshaping_0.3.7    stats4_4.2.2        
## [91] RcppParallel_5.1.7   munsell_0.5.1        viridisLite_0.4.2   
## [94] bslib_0.9.0